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NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
SBIR PHASE H FINAL REPORT 

PROJECT SUMMARY 


CONTRACT NO.: NAS8-39370 


PROJECT TITLE: Radiation from Advanced Solid Rocket Motor Plumes 


PURPOSE: The overall objective of this study was to develop an understanding of solid rocket 

motor (SRM) plumes in sufficient detail to accurately explain the majority of plume radiation test data. 
Improved flowfield and radiation analysis codes were to be developed to accurately and efficiently account 
for all the factors which effect radiation heating from rocket plumes. These codes were to be verified 
by comparing predicted plume behavior with measured test data. Extensive radiation data were provided 
by the NASA/MSFC test programs conducted to support the design of the ASRM. 


RESEARCH COMPLETED: Upon conducting a thorough review of the current state-of-the-art of 

SRM plume flowfield and radiation prediction methodology and the pertainent experimental data base, 
the following analyses were developed for future design use. 

• The NOZZRAD code was developed for preliminary base heating design and A1 2 0 3 particle optical 
property data evaluation using a generalized two-flux solution to the radiative transfer equation. 

• The IDARAD code was developed for rapid evaluation of plume radiation effects using the spherical 
harmonics method of differential approximation to the radiative transfer equation. 

• The FDNS CFD code with fully coupled Euler-Lagrange particle tracking was validated by 
comparison to predictions made with the industry standard RAMP code for SRM nozzle flowfield 
analysis. The FDNS code provides the ability to analyze not only rocket nozzle flow, but also 
axisymmetric and three-dimensional plume flowfields with state-of-the-art CFD methodology. 

• Procedures for conducting meaningful thermo-vision camera studies were developed. 

• The final report on this study provides user’s manuals for the codes developed, and the source codes 
were delivered to NASA for their use. 


RESULTS: The NOZZRAD code was validated for preliminary base heating design use. The FDNS 

code was validated for SRM nozzle analysis. The potential of the IDARAD code was identified for future 
use. New treatments of A1 2 0 3 optical property data base were recommended for making better use of 
existing test data. 
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POTENTIAL APPLICATIONS: The radiation analyses methods developed in this study are useful 

for predcicting and understanding thermal environments of launch vehicles and launch stand facilities, 
of missile infrared missiles, and of decoy designs to defeat heat seeking missiles. The methodology is 
also appropriate to furnace design and the thermal loads produced within gas turbine and diesel engines. 


NAME AND ADDRESS OF CONTRACTOR: 

SECA, INC. 

3313 Bob Wallace Avenue 
Suite 202 

Huntsville, AL 35805 


PRINCIPAL INVESTIGATOR: 

Richard C. Farmer 
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Foreword 


The document presents the results of a phase II SBIR study 
performed by SECA, Inc. to investigate Radiation from Advanced 
Solid Rocket Motor Plumes. The study was performed for NASA 
Marshall Space Flight Center under Contract NAS8-39370. The 
NASA/MSFC technical monitor for the study is Mr. Peter R. 
Sulyma. 
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1.0 INTRODUCTION 


An investigation was conducted to develop an understanding of solid rocket motor (SRM) 
plumes which was sufficiently accurate to explain the majority of plume radiation test data. The 
goal of this study was to produce methodology which can be used as a design tool for predicting 
the radiation environment created by the plumes of a launch vehicle. Considering the vast 
expenditures of many government agencies on experimental studies of SRM plume related 
problems, the successful completion of this investigation offers a significant potential cost 
savings. 

Historically, SRM plume radiation has been underpredicted by existing analytical 
methods. This was attributed to unrealistically low values of the imaginary part of the index of 
refraction of A1 2 0 3 and the neglect of the searchlight effect which redistributes interior motor 
radiation into the plume. Grumman’s shock tube experiments to obtain better A1 2 0 3 optical 
property data (Ref. 1.1) and Remtech’s development of a Monte Carlo code to include the 
searchlight effect (Ref. 1.2) allowed significant improvement in plume radiation predictions. 
However, several additional factors which affect particle size and temperature and the need for 
a more efficient radiation prediction code required more study. This investigation was designed 
to provide these improvements. However, several concurrent researches , namely: (1) the 
continued improvement of optical property data (Ref. 1.3), the development of a new heat 
transfer analysis to account for particle/gas heat exchange (Ref. 1.4), extensive measurements 
of particle size distribution in SRM plumes (Ref. 1.5), and recent access to Russian data on 
plume radiation (Ref. 1.6), greatly influenced the final outcome of this study. One of the more 
significant concurrent studies that helped improve the radiation base heating predictions was the 
development of the improved Cycle 2.0 solid rocket motor flowfield methdology which 
incorporates the results of these recent studies into the heating analysis (Ref. 1.7). 
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This investigation addressed: the evaluation of particle and gas radiation data, the 
development of a computationally efficient radiation analysis, and the prediction of the two-phase 
flowfield. Finally, predictions and comparisons to other methodology of specific test data will 
be presented to evaluate the utility and contribution made by this study. 
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2.0 PLUME RADIATION PROPERTIES 


Radiation properties of the optically active plume species and the method of analyzing 
the radiative transfer process control the accuracy of predicted plume radiation. These subjects 
are described in Sections 2.1, 2.2 and 2.3. 

2.1 Radiation Properties 

SRM plume radiation heating is predominately from A1 2 0 3 particles and is augmented by 
gaseous radiation from combustion products, namely: C0 2 , H 2 0, CO, and HC1. Soot particles 
may also contribute to the radiation, but an accurate measure of such radiation has not yet been 
established. Exhaust products from liquid rocket motors (LRM), which utilize RP-1, H 2 , and 
0 2 as propellants, comprise a subset of these radiating species, hence an analysis which is 
acceptable for SRM’s will also be appropriate for most launch vehicle design purposes. These 
plume constituents radiate in the following manner in the near infrared region of the spectrum 
which controls radiation heating: 

1. A1 2 0 3 and, perhaps, soot particles are large enough that they emit, absorb, and scatter 
radiation in a continuous frequency spectrum. The real and imaginary parts of the index 
of refraction are obtained from experiments, and Mie theory is usually used to convert 
from index of refraction values to absorption and scattering coefficients. 

2. Gaseous radiation in the infrared is non-continuous, non-luminous molecular radiation 
associated mainly with rotational and vibrational energy modes. 

Radiation predictions involve solving the radiative transfer equation (RTE), an energy 
balance on radiation, which requires specification of absorption and scattering coefficients. This 
section addresses the evaluation of such coefficients. 
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2.1.1 Particle Radiation 

Thermal radiation from SRM plumes occurs between 0.5 and 8.0 /xm, as shown by a 
typical SRM spectral signature in Fig. 2.1. The spectral data shown in this figure were taken 
by Sverdrup Technology personnel (Ref. 2. 1) for MNASA-6, a subscale ASRM test motor. The 
spectrometer sees a 6-inch circle at the plume centerline about 24 inches downstream of the 
nozzle exit. The spectrometer was located 340 feet from the nozzle at an elevation of 10 feet. 
The nozzle exhausted upward from an elevation of 17 feet. The radiation peaks in the 1-2 /im 
wavelength region, and absorption by cool combustion gases along the plume boundary and by 
the atmosphere between the plume and the detector cause the dips in the radiance at the band 
centers. A1 2 0 3 optical properties are authoritatively described by Reed, et al (Ref. 2.2) and 
indicate that liquid alumina has an imaginary index of refraction (N 2 ) of: 

N 2 = 4.66E-4 X 1 33 T' 5 exp{-29420/T} 

independent of impurities in the alumina. For solid, crystalline alumina, the state, whether a 
or 7, and impurity levels have a strong effect on N 2 . Particle samples taken from the centerline 
of an SRM plume, designated Rocket 1, and from the edge of the plume, Rocket 2, were 
measured by Grumman (Ref. 2.3). The Rocket 1 particles were found to have less impurities 
than the Rocket 2 particles and to have lower values of N 2 but still higher than pure alumina. 
Rocket 2 particles in the solid phase had an N 2 which was essentially that of the liquid at the 
melt point. Recently, the use of argon as a carrier gas in these shock tube experiments has been 
questioned (Ref. 2.2). When a C0/C0 2 carrier gas was used, lower values of N 2 were measured 
for pure alumina liquid particles. Sufficient experiments to fully qualify this observation with 
respect to particles collected from SRM plumes have not yet been performed. The effect of y 
to a phase transformations will be considered subsequently. For the present, the current 
Grumman shock tube data for alumina particles will be used as the best available data for N 2 
(Ref. 2.2). 


2-2 



'VASA 


SECA-FR-94-18 



2-3 


Fig. 2. 1 Spectra] Radiance from the MNASA6 SRM 
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The real part of the index of refraction, N 1} is insensitive to impurities and to 
temperature. In the spectral interval of interest, it varies slightly with wavelength, and is well 
represented by: 

N, = 1.75 cos{6\} 

where the 6 has units of degrees/micron. 

Optical properties for A1 2 0 3 were obtained from Grumman ’s OPTROCK data (Ref. 2.4). 
Plots of the refractive index N, and the absorptive index N 2 are given in Figs. 2.2 and 2.3. 
These data were transformed to radiation properties by a MIE scattering code (Ref. 2.5) which 
determines absorption and scattering cross-sections as well as phase functions. The OPTROCK 
data includes wavelengths ranging from 0.2 to 25 nm. Notice that between 0.5 and 8 /zm there 
are much less variation in the data. Additional work has been reported (Ref. 2.2) to provide the 
previously mentioned curve-fits of the absorptive index for the liquid A1 2 0 3 optical property data. 
Since A1 2 0 3 radiation is continuum in the region of interest, no problem exists with respect to 
defining a meaningful spectral average for a given wavelength. 

Scattering angle dependence on the diffractive, reflective and absorptive processes are 
determined through the phase function, P, as illustrated in Fig. 2.4. Note that a specific value 
of logj 0 P is computed for each scattering angle COS (0). These figures indicate that as the size 
parameter (X = xD/X) increases, forward scattering increases. The smoothing effect of using 
a particle size distribution is not illustrated by Fig. 2.4, but it is expected. It is known that Mie 
scattering theory is a good approximation for forward scattering for large non-spherical particles. 
Since extinction is dominated by scattering in the forward direction (as seen in these figures), 
then scattering is not very sensitive to particle shape. This justifies using the Mie theory. The 
OPTROCK data for Nj and N 2 were used to create the data tables for a c , albedo, and the 
backscattering fraction in SIRRM (Ref. 2.6). SECA’s MIE code is used to predict a„ and 
the backscattering function. Using both SIRRM and the MIE code results in the radiation 
predictions for a homogeneous slab shown in Fig. 2.5. Data at temperatures of 2300 and 3000 
°K are in the tables; the intermediate temperatures of 2500 and 2800 °K are not. Obviously, 
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differences in the interpolation schemes for the MIE and SIRRM codes exist; methods of 
reducing these differences will be discussed later in this section. 

Several investigators have postulated that rapid quenching of molten A1 2 0 3 particles 
leaves them in a 7 state rather than the more usual a state and that this state exhibits greatly 
different radiative properties. Russian investigators (Ref. 2.7) made a study of this effect with 
both laboratory and flight experiments, but the investigation was limited to small wavelengths, 
such that IR heating could not be evaluated. The Russian investigators experimentally 
demonstrated the importance of alumina impurities and 7 la phase transitions, but offered no 
general purpose models for predicting these effects. The largest wavelength observed was 
1. 1/xm where at most a 20% increase in radiance was observed. PSI investigators (Ref. 2.8) 
observed high emissivities of pure alumina solid particles which were rapidly cooled in a shock 
tube, which substantiates the Russian experiments. However, most of the plume radiance data 
which are difficult to evaluate are for conditions where the particles are molten as they leave the 
nozzle in the plume where the radiance was measured. Sverdrup investigators (Ref. 2.9) 
predicted particle states for a specified phase transition kinetics rate; no optical data were 
presented to support this analysis. Other Sverdrup investigators (Ref. 2.10) sampled rocket 
plumes and determined the ratio of a /7 crystals in the particles and used these data to deduce 
a kinetics expression for the phase transition. As mentioned previously, impurities in the A1 2 0 3 
dominate the solid phase optical properties, and none of the phase transition studies have yet 
addressed the effects of impurities on observed optical properties. 

The data shown for N 2 of solid A1 2 0 3 in Fig. 2.3 suggest that using properties for solid 
particles at the melting temperature is valid at all temperatures; since the room temperature data 
are for pure A1 2 0 3 . This idea has not yet been evaluated with a radiation heating analysis, but 
it should be. 

Although several questions concerning A1 2 0 3 optical property data have not yet been 
satisfactorily answered, available Grumman data are believed to be sufficiently reliable to be 
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used for plume radiation analysis. This is especially true for the purpose of making validation 
predictions for available rocket plume test data. 

Absorption coefficient data for soot are found in the radiation handbook (Ref. 2.9) and 
complex index of refraction values are found in documentation for SIRRM (Ref. 2.6). Other 
data sources for soot are not considered sufficient for consideration. Optical properties for soot 
are reasonably well known; the problem is that the accurate predictions of particle size and 
number density for the soot particles cannot yet be made. 

2.1.2 Gaseous Radiation 

To solve the RTE for gaseous (molecular) radiators, a value of the absorption coefficient 
(jtx) must be provided. For monochromatic radiation this is possible, but too many individual 
lines must be considered to make a practical heat transfer analysis by simply summing the lines. 
Three alternatives have been proposed in the literature: (1) the use of narrow band models (Refs. 
2.6 and 2.11), the use of wide-band models (Refs. 2.12 and 2.13), and the use of a total 
emissivity (Ref. 2.14). It should be noted that linear absorption coefficients (k*) are required 
to solve the RTE; hence, if mass (* pX ) or pressure (* pX ) absorption coefficients are taken from 
correlation equations, they must be converted before they are used. Both the narrow and wide 
band models involve spatial averaging along a line-of-sight (LOS) before spectral averaging over 
a spectral interval can be accomplished. Unless this procedure is reversed, severe restrictions 
on the method of solving the RTE result. Theoretically, the narrow band models involve 
summations over many narrow spectral intervals, whereas the wide band models treat an entire 
band at one time. However, the entire band modeling procedure has not yet been developed to 
the point that a rigorous method of inverting the spatial/spectral integration is available. The 
total emissivity method has only been studied for the C0 2 /H 2 0 system, is completely empirical, 
and is not practical for the several plume species with strong temperature, concentration, and 
pressure gradients in the flowfields of interest. 
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The extensive radiation data developed by MSFC in support of the Saturn program (Ref. 
2. 1 1) provides the narrow band (NB) absorption coefficients for plume gases for many SRM and 
LRM. Inhomogeneous optical paths are treated with the modified Curtis-Godson approximation 
to allow spatial integration along a LOS. Remtech’s Monte Carlo (Ref. 2.15) and GASRAD 
(Ref. 2.16) codes and the SIRRM two- flux and six-flux radiation models use these narrow band 
models for evaluating plume radiation (Ref. 2.17). Integrations along specific LOS are 
performed for spectral intervals of 100 to 400 cm' 1 wave numbers. The narrow spectral intervals 
require lengthy computation times to evaluate radiative heating. Furthermore, local absorption 
coefficients are not provided by the narrow band model which are useful for obtaining solutions 
to the RTE by more general and economical methods. Such absorption coefficient evaluation 
methods could be developed, but the narrow spectral interval integrations would still be required. 

Since the radiation model developed by this study was to be designed to be fast and 
economical to use, the exponential wide band (EWB) model (Ref. 2. 12) was investigated for use 
in solving the RTE. This model was developed to describe radiation from a LOS with constant 
temperature, pressure, and composition for each major band of the optically active species 
present. In general, the test data that the EWB model is based upon were taken at higher 
pressures than the NB model data previously mentioned, but still at pressures much lower than 
typical rocket motor combustion chamber pressure. Correlation parameters for this model are 
given in Table 2.1 for the following thermally important bands for H 2 0, C0 2 , and CO: 

X m o = 1.38, 1.87, 2.7, and 6.3 urn 
Xco 2 = 2.0, 2.7, 4.3, 9.4, and 10.4 /xm 
k co = 2.35 and 4.7 /xm 
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Table 2.1 Wide band model correlation parameters for various gases 



Correlation Parameters 

<*o 

[cm l /(g/m 2 )] 

To 

O ) 0 

[cm 1 ] 


H 2 0 m 

— 3 , t/| — 3652 cm' 1 , i ? 2 

— 1595 cm' 1 , i) 

3 — 3756cm' 1 , gj. 

Ifllllll 


6.3jum 

ij c = 1600 cm 1 
(0,1,0) 

1 

8.6 (T 0 /T) 05 
+0.5 

41.2 

0.09427 

56.4 

2.7/tm 

r) c = 3760 cm' 1 
(0,2,0) 

(1,0,0) 

(0,0,1) 

1 

8.6 (T 0 /T) 05 
+0.5 

0.2 

2.3 

23.4 

1.03219 

60.0 

1.87/xm 

r/ c = 5350 cm' 1 
(0,1,1) 

1 

8.6 (To/T) 1 5 

+ 1.5 

3.0 

0.08169 

43.1 

1.38^im 

T 7 C = 7250 cm' 1 
(1,0,1) 

1 

8.6 (Tq/T) 1 5 
+ 1.5 

2.5 

0.11628 

32.0 

C0 2 m 

— 3, rh — 1351cn 

f = 

666 cm' 1 , i) 3 = 

2396' 1 , g k = (1,2.1) 


10.4/im 

t) c = 960 cm' 1 
(-1,0,1) 

0.8 

1.3 

2.47X10' 9 

0.04017 

13.4 

9 A pm 

ij c = 1060 cm' 1 
(0,-2, 1) 

0.8 

1.3 

2.48X10' 9 

0.11888 

10.1 

4.3/xm 

i] c = 2410 cm' 1 
(0,0,1) 

0.8 

1.3 

110.0 

0.24723 

11.2 

2.7/xm 

7j c = 3660 cm 1 
(1,0,1) 

0.65 

1.3 

4.0 

0.13341 

23.5 

2.0fim 

— 5200 cm' 1 
(2,0,1) 

0.65 

1.3 

0.060 

0.39305 

34.5 

CO m 

- - 2143 cm 1 , g, = 

= 1 


4.7/tm 

jj c = 2143 cm' 1 
(1) 

0.8 

1.1 

20.9 

0.07506 

25.5 

2.35fcm 

i) c = 4260 cm' 1 
(2) 

0.8 

1.0 

0.14 

0.16758 

20.0 



2-12 
























































































SECA-FR-94-18 


Methods to average the EWB model for inhomogeneous path lengths, simply average the 
temperature and species over a finite length (Ref. 2.18). This is not an appropriate method for 
obtaining local values of the absorption coefficient. Modest devised a method for defining local 
coefficients by using an optically thin and a mean beam length to define the two variables: 
absorption coefficient and effective band- width (Ref. 2.19). Even though the accuracy of such 
a method has not yet been established, this method was used by Modest in conjunction with a 
first order ordinary differential approximation. SECA used the same EWB model in the 
ordinary and improved differential approximation methods which are described in the next 
Section. 

Ultimately, both particle and gaseous radiation must be described. Figure 2.6 shows a 
LOS calculation for a particle gas mixture which was made with the narrow band models in 
SIRRM (Ref. 2.6). This example was taken from a slice out of a SRM plume. The radiation 
looks like an averaged absorption coefficient for the A1 2 0 3 particle gas mixture was used. 
However, gas/particle radiation interaction can look quite different, depending on specific 
concentration and temperature variations along the LOS. 

The narrow band models, such as those used in the JANNAF standard plume radiation 
code, SIRRM, are the most accurate molecular radiation models currently available. However, 
the spectral integration required to use such models is too computationally intensive for 
economical use and probably prohibitive for practical use in a fully coupled solution. Therefore, 
an evaluation of the EWB model was made in this investigation. 

2.2 Solution of the RTE by Using Differential Approximations 

The radiation analysis developed for this study is uncoupled, in the sense that the 
flowfield is first calculated then the radiation resulting from the predicted temperature 
distribution is calculated. Such a treatment assumes that the energy lost by radiation is small 
compared to the energy in the flowfield. The plume’s of SRM’s emit, scatter, and absorb, so 
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Fig. 2.6 Gas/Particle Radiation Interaction 
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that a solution to the general radiative transfer equation (RTE) is needed for the intensity, I x . 
Hence the intensity obtained from the RTE given below must be integrated in wavelength to 
obtain heating rates. 

• VX A {f, 0} = K A J Ai> {T} - (k a + a x )I x {f,Cl} 

a (1) 

Jp{?, fi'-fi} I x {f, £J'} dQ' 

4 71 

The RTE may be formally integrated over space to obtain: 

!{?, £5} = J{f\ £5} e' T 


♦/ 







J{f, £}} P {r, fi'-fijdfi' 


dt 


( 2 ) 


The subscript X is suppressed in this and subsequent equations for clarity, but it must be 
remembered that the intensities are appropriately averaged monochromatic values. Direct 
numerical solutions to the integrated RTE have been reported by Tan (Ref. 2.20) and Tan and 
Howell (Ref. 2.21). These methods are quite interesting, but such solutions have not been 
performed for problems as complex as those found in rocket plumes. To provide efficient 
solutions to the RTE, a variation of a differential approximation and a generalized two-flux 
model were chosen for further development in this study. Details of these methods are described 
in the remainder of this Section. Of course, the Monte-Carlo method (Ref. 2.15) could have 
been used for solving the RTE, but this method was not deemed to be sufficiently fast for 
exploratory studies or for routine use. 


2.2.1 The ODA Method 


The method of spherical harmonics (a type of differential approximation) was developed 
for solving the RTE. The method of spherical harmonics is implemented by expressing the 
phase function and intensity as Legendre polynomials as shown in Chart 2.1. The result of these 
transformations is to replace the integro-partial differential equation with partial differential 


2-15 



SOLUTIONS TO THE RTE USING SPHERICAL HARMONICS 


SECA-FR-94-18 



2-16 


Limitation: For N=1 or 3, accurate for only optically thick media. 
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equations (PDE). The P, approximation or the ordinary differential approximation (ODA) is 
shown along with the one PDE which must be solved, but this method is accurate only for 
optically thick media. Higher order P N approximations may be used, but the accuracy increases 
slowly with the order of the solution whereas the number of partial differential equations 
required for a solution increases dramatically. For the P 3 approximation 4 simultaneous PDE 
must be solved. Modest (Ref. 2. 19) has suggested the Modified Differential Approximation 
(MDA) and the Improved Differential Approximation (IDA) as shown in Charts 2.2 and 2.3, 
respectively. The IDA is simpler to use, and it has the advantage of starting from the P, 
solution. Hence, if the media is optically thick, i.e. inside the motor, only this part of the 
analysis is required. Having the P, solution, the J, G, and (f terms shown in Chart 2.3 must 
be evaluated for each point where a solution is desired, i.e. at each grid point in the flowfield. 


To obtain the Pj solution, the RTE will be written explicitly in terms of cylindrical 
coordinates (Ref. 2.22) as 


cos sin0-^- + sin (<j)~4> ) sin6 + cos Q-^- 

dr i o<|> r dZ 

+ J v (r,<|> r ,Z;0,<|>) = (!-«) I 


Z;0,< t>) 


2n 1 


f f T v (r, 4> r , Z; 0' , <j)' ) P (0, <J>;0' , <J>' ) d cos 0 ' d<{>' 


O -1 


(3) 


To transfer Eq. 3 to the P, -approximation, the phase function is defined by: 

(d) yT id') (4) 

1*0 m=~l 


Likewise, I, is expanded in spherical harmonics. Next, the first moment of the reformulated 
RTE is made, and after the same manipulation the P,-approximation in cylindrical coordinates 
results in a Helmholtz equation for I 0 „: 
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THE MODIFIED DIFFERENTIAL APPROXIMATION 

The MDA generalizes P, solution to treat arbitrary optical thickness. 
Let surface & media contribute to intensity: 
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+ j!£°i _ a 2 j 

r ar^ ar J r 2 a<j > 2 az 2 


- A 2 r fcv {r} 


( 5 ) 


For an axisymmetric system, the <f> term is zero. 
A 2 = 37 0 7 t aj 

7o = 1 - up,, 


7 , = 1 - cop/3 


(for isotropic scattering, p 0 = 1, p, = 0) 

I*, is the first moment of the scattering integral defined by: 

2 k it 

W/'. sin0d0d<J> (6) 

O O 

The radiation intensity in Eq. 6 is a function of position and direction. However, the evaluation 
of the Helmholtz equation above for the P, -approximation eliminates the need to resolve the 
complicated angle dependencies within the integral. The azimuthal angle <£ was eliminated from 
the original RTE by integration over all directions. In other words, when the first moment of 
the RTE is taken to derive the P,-approximations, the <f> term becomes: 


1 dl o V 

r 3(f> r 


2w * 

// sin (<|>-<l) r ) sin 2 0 dOcftj) 

O O 


1 


r 



= 0 


(V) 


Equation 5 has been written in a discretized form suitable for computer simulation. The 
procedure for solving Eq. 4 proceeds from a successive-line overrelaxation (SLOR) method. 
Rewriting Eq. 5 with a source term defined by R IOU , 
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v^J ov = A 2 = R Im (8) 

V c 2 is the cylindrical operator corresponding to the differential terms on the left-hand side of Eq. 
5. 

Terms in Eq. 8 have been normalized with respect to reference values and will not be 
discussed further. 

( symmetry ) r=r, 

(inlet) Z=Z 1 (~) ^ 

(exit) Z=Z 2 (+) 

where h = 1.5a ext 7,. At the plume boundary there is no incoming radiation. 

Once the P, solution is obtained, the IDA analysis can be accomplished. Since this 
radiation analysis is to be uncoupled, provision for using either SPF/2 and FDNS flowfield 
solutions will be provided as input for the plume radiation analysis. 


Boundary conditions for Eq. 5 are: 


dl. 


OV _ 


di 


dr 


= 0 


3 #* * 0 


2.2.2 The IDA Method 


Since the ODA is not valid for optically thin regions, it was extended to accommodate 
all optical thicknesses through the IDA. The logic which allows the ODA to be transformed 
to the IDA lies in the linear approximation of the radiative source term s*, which is itself a 
function of the ODA results. Also, the radiation intensity at any point (i, j) in the medium is 
written as the sum of wall and medium contributions. Transmissivity of wall radiosities to the 
medium point are obtained from extinction coefficients that are evaluated from a Mie scattering 
code. The medium contribution includes an adjusted source term which is not evaluated at the 
medium point (i, j), but at an adjusted point. Instead of evaluating the derivative of s* at (i, j), 
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it is easier to evaluate S at a point r 0 away from (i, j). This is accomplished by the following 
procedure. 

A Taylor series representation for s can be written at a point (i, j) or T = 7 W + s^ in 
the medium by assuming that s varies linearly from T to 7 W in the direction (-§) by: 

+ s w §, 8) = s*(f, §) - (t s -t ( 10 ) 

r, is the optical depth along a path length from (i, j) to 7 W , the star refers to values based upon 
ODA results, and subscript w corresponds to the wall. S w is the magnitude of the vector § from 
the wall to r w . The radiative intensity in the medium is determined by straight-forward 
integration of: 


= /«*(*„ + sj, 8) dr' 

o 

( 11 ) 

= s*(?,§) (2-e~ T *) - — [i - (1 + x s ) e _t *] 

^ * a 


To eliminate the derivative in Eq. 1 1 , the method is to take the assumed linearity of s in reverse 
order, transforming an integral to an unknown Taylor series. The procedure is to determine the 
Taylor series, call it L, that produces the RHS of Eq. 11: 

L = RHS ( 2 ) (12) 


Rearranging terms in Eq. 12 leads to 


L 

l-e Xa 


8 ) 



1-e' 


\ 


/ 


ds * 
d's 


(*, 6 ) 


Comparison on Eqs. 10 and 13 gives the following correspondences: 


( 13 ) 
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-T.\ 


1-e' 




S*(f w + s w §, §) 


In this way, (7, - 7, w ) is replaced by 7 0 , where 


t =1 - 


l-e" 


( 14 ) 


Then (s^) is replaced by (-s 0 s) since the series is taken in reverse order, so that the radiative 
source term becomes s*(r w -s 0 §,§), and s 0 is backed out from 


— o 

*0 = / o ext {Z-s"§)ds" 


( 15 ) 


Finally, from 

- s-a.sj.g) 

1-e • 

we obtain the result 


J_ = 


s 

L = fs*(f w + s w 8,6)e’ lx ‘- x ' ) dz' s 


( 16 ) 


* s*(f w -s 0 §,§) (i-e“ T *) 
s 0 is the physical distance across the LOS into the plume, and 


S*{r w -s 0 s, s} = (l-oj)I* b {r w -s 0 s } + (co/47r)[G*{r w -s 0 s} + A,q*{r w -s 0 s}* s] (17) 
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Eq. 16 replaces Eq. 10 and does not require an evaluation of a derivative, but the value for s Q 
must be found along each line-of-sight according to Eq. 15. The radiative source contribution 
at (i, j) is transformed to evaluating a radiative source a distance s 0 from point (i, j) in the 
direction toward the wall. 

The IDA modules were coded for accepting FDNS or SPF/2 flowfields as input. Briefly, 
the IDA modules are set up to allow the user to choose the number of wall points deemed 
important for radiosity computations; a reduction in CPU requirements will result if fewer wall 
points are chosen. A view factor code then evaluates view factors between wall segments, after 
which the optical distance between wall points and between medium and wall points are 
determined. Adjusted source terms, wall radiosities, and wall and medium contributions to the 
incident radiation are then evaluated. 

Analyses for gaseous radiation based upon the exponential wide band (EWB) model for 
the absorption coefficients of H 2 0 and C0 2 was included in the ODA and IDA codes. The 
rocket fuels of interest create aluminum oxide particles and soot in the combustion process and 
other gases, primarily CO and HC1. Using the FDNS-EL code, two-phase flowfield predictions 
can be post-processed through a SIRRM map module that determines the particle density, gas 
phase species concentrations, pressure, and temperature as a function of position. Combining 
the Mie code with the optical property data for A1 2 0 3 to get absorption and scattering 
coefficients, along with the SIRRM map, particle radiation computations can be made. This 
procedure has been implemented for the ODA where the scattering integral is now accounted 
for through the albedo and an isotropically scattering phase function. Multiple solutions of the 
RTE equation are used to predict monochromatic or band averaged intensities; these intensities 
are then summed to obtain total radiative transfer. 

A similar procedure was used for the IDA. The only remaining step for the IDA 
gas/particle radiation method was the development of an optical path procedure between two 
points so that surface radiosity and source effects can be determined. Since the ODA method 
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is valid only in optically dense regions and can be substantially in error when surface emission 
dominates over medium emission, it was extended to include arbitrary optical densities through 
the IDA method. This is especially necessary for plume regions where the optical thickness is 
relatively small. The main concept in the IDA is to bring in the wall or geometric influence to 
the radiation calculation at all flowfield points. In the ODA methodology, the Legendre 
Polynomial Series expansion for the radiation intensity and the phase function were sufficient 
to account for an optically dense region. But in optically thin cases, the influence of radiation 
from one point to another requires extending the optical path farther out. The IDA procedure 
accomplishes this by incorporating the wall effects into the formulation. It is therefore important 
to be able to choose, for each flowfield node, a sufficient number of lines-of-sight to the wall 
surfaces to properly account for the radiosity effects. To simplify the logic, the radiation 
intensity is split into two terms, one for the wall and one for the medium. The three steps to 
an IDA solution are: 1) ODA solution for the flowfield node source term, 2) surface integrals 
for wall radiosities, and 3) surface integrals for flowfield node incident radiation. 

An important part of the surface integral, the geometric component, requires the 
evaluation of the view factors between two wall points and of the solid angles between a wall 
and a flowfield node. An existing code, RAVFAC (Ref. 2.23) was incorporated into the IDA 
solver, along with a preprocessor code that initializes RAVFAC with surface data for a nozzle 
configuration. The coordinates of the nozzle wall are obtained from a grid file, and then the line 
connecting two neighboring wall boundaries are described as either a circular disk (inlet), a 
cylinder (combustor wall) or a cone (converging and diverging nozzle sections). Other input 
variables describe the local coordinate system of each surface shape to allow the evaluation of 
unit normals and the determination of whether a surface is shaded by another surface. The 
accuracy in describing view factors between wall surfaces can be chosen relative to the number 
of angular planes desired. This can range from a single 360-degree circumferential surface, 
resulting in circular band or hoop shapes, to any portion of 360 degrees, which results in a much 
larger number of wall-to-wall combinations. In addition, code was provided to allow the 
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calculation of the solid angles between wall points and flowfield points, which is needed for the 
evaluation of the incident radiation. 

Another important aspect of the IDA requires the evaluation of source terms and optical 
distances along various lines-of-sight. As mentioned previously, the flowfield influence to the 
IDA equations includes the source term S* from the ODA as part of its solution. Instead of 
using S* at a particular flowfield point, however, an adjusted S* is required. A module that 
determines the optical length along a chosen line-of-sight was added to the IDA code, along with 
the logic that determines the adjusted S\ Lines-of-sight are chosen to extend from the 
axisymmetric plane, where the flowfield solution is known, to other circumferential planes. The 
line-of-sight module evaluates the product of the extinction coefficient and a differential distance, 
for each increment along a line-of-sight from a flowfield node to a wall node, and sums these 
values for a total optical depth r,. An adjusted source location S 0 along the line-of-sight is 
backed out from t 0 = fn (r,). Once the value of S 0 is determined, its coordinates are extended 
back to the axisymmetric plane, from which the new S* can be interpolated. The surface 
integrals for radiosity, J w , and incident radiation, G, then follow from the variables described 
above. An inversion routine for J w that is efficient in inverting a matrix with non-zero entries 
in almost every location (the entries have a zero value where the view factors between surface 
nodes are zero) was also included. Specifically, the equation for the wall radiosity can be 
written as: 



£ 

ii 

✓a 

(18) 

where the elements of matrix A are 



a u 

= l-CL-c*) e~ %ii F n 

(19) 

a u 

= -(!-€_£ ) e' T « F i:i 

(20) 


and the right-hand-side is 
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R i ~ e i n I bi + U-e'^) F iJ (2i) 

J 

The Fjj is the view factor from surface node i to surface node j. This matrix differs from the 
usual equation for J w by the addition of the adjusted ODA source term S*. 


The discussion to this point has focused on the radiation solution within a region enclosed 
by solid walls, an inlet and an exit. The method was extended to plume regions where radiation 
heating to a rocket base region must be determined. To accomplish this task appropriate 
boundary conditions must be specified on the boundaries for the initial ODA solution. In 
particular, if the plume region is solved independently from the nozzle flowfield (after the nozzle 
solution is obtained), the nozzle exit can be prescribed with the following condition: 


J = 


71 


( 22 ) 


with the percent of radiation crossing the exit plane specified by a (Ref. 2.24), and with an 
emissivity specific to the gas and particle mixture at each exit point. For the other boundaries, 
a constant radiation intensity can be specified as above. Lines-of-sight from the rocket base 
region to any point within the plume can be accomplished as was done within the nozzle. 


Initially, the ODA equations were solved for a specified wavelength, where the 
wavelengths corresponded to those of the gas or particles. The ODA equations were then solved 
for each wavelength without regard for gas/particle overlap. Logic was added to extend this 
procedure by taking into account the fact that gas band radiation is modified by the presence of 
the overlapping particle bands. Specifically, the gas band absorption coefficient was augmented 
by the extinction coefficient of the particle background (Ref. 2.19). A set of subroutines was 
developed to separate the particle-only bands or windows from the gas/particle overlap bands. 
The option to run a single wavelength for purposes of getting an initial solution is available as 
well, along with the option of running a monochromatic solution with a mean absorption 
coefficient for H 2 0 and C0 2 . 
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A listing of the types of options available in IDA is given below in Table 2.2, followed 
by a brief description of the subroutines and a flow chart of the code, which are shown in Table 

2.3 and Chart 2.4, respectively. The BLKDAT subroutine contains all of the radiation 
parameter options necessary to make an ODA or IDA computation. The 9 major options shown 
in Table 2.2 must be chosen to initiate the radiation computation. 

2.2.3 Description and Use of IDARAD Code 

The overall architecture of the IDARAD radiation code evolved as a result of using a 
version of the computational fluid dynamics (CFD) code, FDNS (Ref. 2.25), that was developed 
under a previous NASA phase I SBIR study (Ref. 2.26) to examine the coupling of the radiation 
and fluid mechanics that can occur in high temperature, high pressure rocket exhausts. Input 
and use of the resultant radiation code is somewhat more cumbersome than would otherwise 
have resulted if the code had been written from scratch. However, a large advantage of using 
the FDNS code is the ability to incorporate the ODA and IDA methodologies in a coupled CFD 
code so that at some future date, a fully coupled radiation/fluids code could be more easily 
developed. 

A flow chart of the program elements that make up the IDARAD radiation code along 
with the data files and transfer of information between the various codes is shown in Chart 2.5. 
The code generates the grid (fort. 12), and flowfield property (fort. 13) and mapped particle-gas 
property (fort.61) files. The RADO program generates initial guesses for the radiant intensities 
of the individual gas and particle bands (fort. 15 and fort. 16) at each of the grid locations. 
Following the execution of these codes, the IDARAD code is executed to calculate the heat 
fluxes at the boundaries. Details of preparing the user specified input files and preparing the 
codes for execution is found in Section 2.2.3. 1. 
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Table 2.2 IDA Code Options 


1. 

Radiation wavelength input: 

• gas species type 

• particle type (presently A1 2 0 3 ) 

• values of wavelengths of interest for gas and particle species 

2. 

IDA switch 

3. 

View factor and solid angle switch 

n 

Wall and inlet plane gas emissivity 

5. 

Spectral or mean absorption coefficient switch for H 2 0/C0 2 mixture 

6. 

Absorption coefficient and band width model: 

• exponential wide band 

• picket fence 

• box model 

• block model 

1 

Mean beam length for cylinder, function of: 

• average distance between neighboring grid points 

• diameter at given axial location 

• radius and length 

• radius and optical depth at band head 

8. 

Boundary conditions at inlets and wall: 

• spectral diffuse emitting and reflecting wall 

• diffuse emitting and reflecting BC with pseudo-black inlet 

9. 

Gas/particle radiation procedure: 

• gas only 

• particle only 

• gas + non-overlapping particle 

• gas/particle overlap 


















Table 2.3 Description of IDA Code Subroutines 
(in alphabetical order) 
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ASTAR 

BLKDAT 

CELLU 

ECOEFG 

ECOEFP 

EXTINC 

GMEDUM 

GWALL 

HTRAD 

IDA 

IDABC 

INRAD 

INTRP1 

INTRP2 

INTRP3 

JWALL 

LININR 

LININT 

LOGINT 

LSIGHT 

MONCHR 

PDA 

ODABC 


Dimensionless band absorption for H 2 0 and C0 2 
Radiation parameter options 

Determines cell in which adjusted IDA source term lies 

Extinction coefficient for H 2 0 and C0 2 

Extinction coefficient for A1 2 0 3 

Driver for extinction coefficient 

Driver for medium surface integrals 

Driver for wall surface integrals 

Radiative heating 

Driver for IDA method 

Determines incident radiation on radiating boundaries 

IDA value for incident radiation and radiative heat flux 

Interpolation of incident radiation from IDA grid to total flowfield grid 

Interpolation of SIRRM properties along a line-of-sight 

Log 10 interpolation of adjusted IDA source term within a cell 
Sets up matrix elements for radiosity 

Linear interpolation of refractive index as function of te mperature 

Double linear interpolation of particle size and temperature for the 
absorption coefficient 

Double log 1 0 interpolation of particle size and temperature for the 
scattering coefficient 

Determines coordinates on a line-of-sight through a 2D axisymmetric 
geometry 

Monochromatic values for absorption and extinction coefficients, albedo 
and Planck function 

Driver for ODA method 

Boundary condition module for ODA incident radiation 
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Table 2.3 Description of IDA Code Subroutines (Continued) 
(in alphabetical order) 


OPDIST 

Driver for optical depth along a line-of-sight 

OVRLAP 

Interpolation for particle absorption and scattering coefficients 
within a gas band 

PATH 

Path length used for gas extinction coefficient 

PBAND 

Particle band widths 

PLANKF 

Planck blackbody function 

RAD0 

Driver for radiation code 

RADIN0 

Driver for radiative transfer equation 

RAVFACSA 

Solid angle routine 

RAVFACVF 

View factor routine 

RDOSTY 

Radiosity matrix inversion routine 

REFRAC 

Reads in particle refractive index (used in evaluation of Planck 
function) 

REFRIN 

Driver for refractive index interpolation 

RSIRRM 

Reads in SIRRM map 

SIGACL 

Reads in particle absorption coefficient data file 

SIGSCL 

Reads in particle scattering coefficient data file 

SIGAV 

Driver for particle absorption and scattering coefficient 
interpolation 

SRCIDA 

Driver for IDA (source terms adjusted for transmissivity and 
absorptivity) 

SUMMA1 and 
SUMMA2 

Wide band model summation - function of vibrational quantum 
number 

TAUDR 

Determines optical depth at single position on a line-of-sight 

TRDIAG 

Tridiagonal solver for ODA 

VFINIT 

Initializes RAVFAC input file for a nozzle 

WIDEBM 

Wide band model for H 2 0 and C0 2 
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Record 1: 
Record(s) 2: 

Record (s) 3: 

Record(s) 4: 

Record(s) 5: 

Notes: 


Table 2.4 Description of fort.12 Grid File 


(Free Format) 

IZON Number of grid zones 
(Free Format) 

IZT(IZ) Number of i stations (x direction) for zone IZ 

JZT(IZ) Number of j points (y direction) at each station for zone IZ 

KZT(IZ) Number of k points (z direction) at each station for zone IZ 

Format: 5(1PE16.8) 

These records input the axial(x) locations of each grid point. 

(X(I,J,K), 1 = 1, IZT(IZ)), J=l, JZT(IZ)), K=l, KZT(IZ)) 

The X values are non dimensional values based on XREF input in fort. 1 1 file. 

Format: 5(1PE16.8) 

These records input the radial (Y) location of each grid point. 

(Y(I,J,K), 1=1, IZT(IZ)), J = l, JZT(IZ)), K=l, KZT(IZ)) 

The Y values are non-dimensional values based on XREF input in fort. 1 1 file. 

Format: 5(1PE16.8) 

These records input the Z direction location of each grid point. Set = 1.0 if 
asymmetric case is being run. 

(Z(I,J,K), 1 = 1, IZT(IZ)), J = l, JZT(IZ)), K=l, KZT(IZ)) 

The Z values are non-dimensional values based on XREF input in fort. 1 1 file. 

If multi-zone cases are being input, input all records 2 first followed by records 
3,4, and 5 for each zone. The radiation code is limited to one zone while FDNS 
can be run using several zones. 
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Table 2.5 Description of fort. 13 Flowfield File for Radiation Code 
Record 1: Format: 815 


Variable 

Column 

Value 

INSO(l) 

5 

1 

INSO(4) 

10 

1 

INSO(5) 

15 

1 

INSOFM 

20 

1 

NGAS 

24-25 

12 


The variables input on Record 1 control the input for the flowfield properties that are 
input to the radiation code. The radiation code assumes that the gas specie set consists of 12 
species in the order specified on the example fort. 11 file shown in Table 2.6. 


Record(s) 3: Format (5(1PE16.8) Gas Density 

These records input the gas density at each grid location 
(((DEN(I,J,K), 1 = 1, IZT(IZ)), J = l, JZT(IZ)), K=l, KZT(IZ)) 

Densities are input in non-dimensional values based on DENREF input in fort. 1 1 
file. 


Record (s) 4: Format (5(1PE26.8)) 

These records input the component of velocity in the i (x) direction at each grid 
location 

(((U(I,J,K), 1 = 1, IZT(IZ)), J = l, JZT(IZ)), K = l, KZT(IZ)) 

Velocity is input in non-dimension values based on UREF input in fort. 11. 

Record(s) 5: Format 5(1PE16.8) 

These records input the component of velocity in the j (y) direction at each grid 
location 

(((V(I,J,K), 1=1, IZT(IZ)),J=1, JZT(IZ)), K=l, KZT(IZ)) 

Velocity is input in non-dimensional values based on UREF input in fort. 1 1 file. 

Record(s) 6: Format 5(1PE16.8) 

These records input the component of velocity in the k(z) direction at each grid 
location 

(((W(I,J,K), 1=1, IZT(IZ), J=l, JZT(IZ)), K=l, KZT(IZ)) 

Velocity is input in non-dimensional values based on UREF input in fort. 1 1 file 
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Table 2.5 Description of fort. 13 Flowfield File for Radiation Code (Continued) 


Record(s) 7: Format 5(1PE16.8) 

These records input the static pressure at each grid location 
(((P(I,J,K), 1 = 1, IZT(IZ), J=l, JZT(IZ)), K=l, KZT(IZ)) 

Pressure is input in non-dimensional values based on the product 
DENREF*UREF**2 input in fort. 1 1 file 

Record(s) 8: Format 5(1PE16. 8) 

These records input the temperature at each grid location 
(((TM(I,J,K), 1=1, IZT(IZ), J=l, JZT(IZ)), K=l, KZT(IZ)) 

Temperature is input in non-dimensional values based on the reference 
temperature-TREF input in fort. 1 1 file. 

Record(s) 9: Format 5(1PE16.8) 

The records input the turbulent kinetic energy at each grid point. Record 9 is 
input only if INSO(5) = 1 

(((DK(I,J,K),I = 1, IZT(IZ)), J=l, JZT(IZ)), K=l, KZT(IZ)) 

Turbulent kinetic energy is input in non-dimensional values based on the 
reference velocity squared (UREF**2) input in fort. 11 file. 

Record(s) 10: Format 5(1 PEI 6. 8) 

These records input the turbulent dissipation at each grid location, Record 10 is 
input only if INSO(5) = 1 . 

(((DE(I,J,K), 1 = 1, IZT(IZ)), J = l, JZT(IZ)), K=l, KZT(JZ)) 

Turbulent dissipation is input in non-dimensional values based on the product 
UREF**3/XREF input in fort. 11 file. 

Record(s) 11: Format 5(1PE16.8) 

These records input the Mach number at each grid location 
(AP(I,J,K), 1=1, IZT(IZ)), J = l, JZT(IZ)), K = l, KZT(IZ)) 

Mach number is input in non-dimensional values based on AMC input in fort. 1 1 . 

Record(s) 13: Format 5(1PE16.8) 

Records 13 are input only if NGAS>0 from fort. 11 file. Records 13 input the 
mass fractions for each gas specie in the same order the species are input in 
fort. 11. 

(((FM(I,J,K,L), 1=1, IZT(IZ)), J=l, JZT(IZ)), K=l, KZT(IZ)), L=l, NGAS) 
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Table 2.5 Description of fort. 13 Flowfield File for Radiation Code (Continued) 

Record(s) 14:* Format 5(1PE16.8) 

These records input particle temperature for each grid location in dimensional 
values (Deg K) 

(((TMPP(I,J,K), 1 = 1, IZT(IZ)), J=l, JZT(IZ)), K=l, KZT(IZ)) 

Record(s) 15:* Format 5(1PE16.8) 

These records input particle number density at each flowfield point in dimensional 
values (#/CM**3) 

(((DNPP(I,J,K), 1 = 1, IZT(IZ)), J=l, JZT(IZ)), K=l, KZT(IZ)) 

Note:* Records 14 and 15 are input only for uncoupled radiation cases (IFL13 > 0). For 
each particle group (size) records 14 and 15 are input followed by records 14 and 
15 for each size group until all sizes have been input. The uncoupled radiation 
code assumes that 5 particle groups are always used. 
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The starting point for the radiation analysis is a precomputed flowfield file in a standard 
JANNAF SIRRM (Ref. 2.16) format. The SIRRM map describes the spatial variation of gas 
pressure, temperature and species as well as particle temperatures and number densities. 

The basic FDNS code requires 3 input files. Unit fort. 11 is the input data file that 
controls the operation of the FDNS or radiation code. Section 2.2.3. 1 describes the generation 
of, and variables input into this file. The other two files are the grid file (fort. 12) and the 
flowfield file (fort. 13). The format of these two files are shown in Tables 2.4 and 2.5, 
respectively. 

The FDNS code uses grid systems that follow the right hand rule for the i(X), j(Y) and 
k(Z) line orientations. The basic FDNS code can handle two-dimensional (axisymmetric) or 
three-dimensional grids having several zones. However, the radiation code developed under this 
contract is limited to axisymmetric single zone grids. The grid must be represented by i axial 
stations with j radial points at each station. Z values must be set to one. It is not necessary to 
have a large number of grid points. A typical mesh for a sea level plume (out to 10 exit radii) 
would be 20-30 i stations with 20 j points per station. Results for radiation at the boundary of 
this size grid versus a grid having 30 times as many grid points are basically the same. Run 
times for an IDA solution for a 400 grid problem of this type are on the order of 6 minutes of 
CPU time on an IBM 320 RISC system. 

The flowfield file (fort. 13) described in Table 2.5 provides the gaseous and particulate 
properties at each grid point. All variables except the particle temperatures and number densities 
are input non-dimensionally based on the reference values prescribed in fort. 1 1 (Table 2.6). The 
only variables that are actually used by the radiation code are gas pressure, gas temperature, 
specie mass fractions, particle temperature, and number density. All other variables can be set 
to 1.0. 


The user must also input a SIRRM. MAP file corresponding to the same grid as the 
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fort. 12 and fort. 13 files. This is really redundant, but due to how the code was developed and 
coded, makes it necessary to input this file. SIRRM maps from typical flowfield codes such as 
SPF/2 (Ref. 2.27) have non-uniform (not equal number) of radial points as well as having too 
many points. A preprocessing code called GRID has been written that uses the SIRRM flowfield 
the user supplies to generate the fort. 12, fort. 13 and SIRRM flowfield (fort.61) that is input to 
the radiation code. 

The GRID code allows the user to specify how many radial (j) grid points to use for the 
grid and how many axial stations to skip between i stations for the grid. Thus, a number of 
SIRRM map points can be eliminated when generating the grid (fort. 12), flowfield (fort. 13) and 
map (fort.61) files. This is a interactive code that requires the user to have the input SIRRM 
map file named as SIRRM.MAP. The user must respond to inquiries by the GRID code as to 
how many radial points are desired (as well as whether the points should be evenly spaced or 
compressed toward the outer part of the flowfield), the station at which to begin the grid, the 
station to end the grid and reference values of length, density, velocity and temperature. These 
values are used to non-dimensionalize values of the grid and flowfield files and must be the same 
that are specified on fort. 11. The GRID code assumes that there are 12 species in the same 
order as is specified on the sample fort. 11 file shown in Table 2.6 (i.e. H20,02,H2,0, H, OH, 
CO, C02, Cl, C12, HC1 and N2). If any other set of gas species are input on fort. 11, the data 
statement for the specie names in GRID.f must be changed. The source code for GRID and the 
SIRRM map for the MNASA48 ASRM contoured nozzle Cycle 2.0 plume are contained on the 
MS-DOS disk, DISC. 

2.2.3. 1 Preparation of Input Files, Subroutines and Steps Necessary to Run ID ARAD 

Each time IDARAD is executed for a new case, several steps must be performed prior 
to execution. In addition to preparing the grid (fort. 12), flowfield (fort. 13) and SIRRM map 
files (fort.61), which was described in Section 2.2.3, the following steps must be performed: 
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(1) Generate the input data file fort. 11. A description of the fort. 1 1 file is contained 
in Table 2.6. This table describes the fort. 1 1 input for both ID ARAD as well as 
the FDNS code, which was used to produce the FDNS results presented in 
Section 3 of this report. FORTRAN unit 11 (fort. 11) inputs to FDNSEL and 
ID ARAD are slightly different and are noted in Table 2.6. Further explanation 
of these inputs can be found in Ref. 2.25. A sample input file for radiation 
predictions inside a solid rocket motor using the IDA method can be found in 
Table 2.7. 

(2) The second file that must be generated is the NOZZRAD.INP file. This file is 
necessary only if running the IDA model. Table 2.8 presents a discussion of the 
input variables contained in NOZZRAD.INP. Table 2.9 presents a listing of the 
NOZZRAD.INP file that corresponds to the fort. 11 file shown in Table 2.7. 

(3) Modify the radOl include file for both the IDARAD and RADO programs. Table 
2. 10 describes the variables that are contained in radOl . These variables set array 
sizes inside the codes. 

(4) Modify the BLKDAT_U.f (IDARAD) and BLKDAT_0.f(RAD0) subroutines to 
set the proper input parameters that are specified by these routines. Table 2.11 
provides a description of the variables contained in the BLKDAT_*.f files. 

(5) Recompile IDARAD and RADO if either BLKDAT_*.f or radOl were changed 
from previous runs. If radOl was changed, recompile the entire program since 
the radOl file is an ’include’ file called by numerous subroutines. 
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Table 2.6 Description of FDNS and ID ARAD 
(fort. 11) Input Data File 

Record Group #1: Gives the case title and identifies whether the problem is 2-D or 3-D. 

Format: 

IDIM, (put title of the problem here — maximum 60 characters) 

< (one data line) 

Definition: 

IDIM = 2 for 2-dimensional flow problems 
= 3 for 3-dimensional flow problems 

Record Group #2: Specifies zonal information and number of flow and wall boundaries. 


Format: 

IZON, IZFACE, IBND, ID, ISNGL (FDNS) 

IZON, IBND, ID, IRAD, IDRW, IFL13 (IDARAD) 

< (one ii ne ) 

Definition: 

IZON number of zones or mesh blocks 

IZFACE number of patched interfaces 


IBND number of flow boundaries (e.g. inlet, outlet or symmetry 

planes) 

ID number of wall elements (blocks) 

ISNGL number of singularity lines 

IRAD radiation control parameter 

0: No radiation: >0: Radiation 

1: Gas rad only 

2: Particle rad only 

3: Gas and particle radiation separately (with no 

gas/ solid overlap regions) 

4: Gas and particle radiation with overlap (with at 

least one gas/solid overlap region). Treats the 
overlapping particle band with the same band width 
as the overlapping gas band width 

a. Does not account for overlap of adjacent gas 
band and points 

b. Does not account for gas overlap at particle 
band endpoints 

IDRW Number of radiating boundaries in axisymmetric plane (for 

IDA method) (eg, can be solid wall or inlet) 

IFL13 0 Do not input particle properties on Unit 13 file 

1 Input particle properties on Unit 13 file 

* Each card group has a header card of whether the record is used or not. See Table 2.7 for examples for IDARAD 
and Table 3.1 for FDNS. 
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Record Group #3: Specifies zonal grid size and zonal rotational/ translational speeds 

Format: 

IZT, JZT, KZT, CBGX, CBGY, CBGZ, CBVX, CBVY, CBVZ 

< (IZ = l,IZON)(FDNS) 

IZT, JZT, KZT (IDARAD) 

Definition: 

IZT(IZ) I-max in zone IZ 

JZT(IZ) J-max in zone IZ 

KZT(IZ) K-max in zone IZ 

CBGX(IZ) rotational speed (RO x /U rcf ) of zone IZ about X-axis 

CBGY(IZ) rotational speed (RO y /U ref ) of zone IZ about Y-axis 

CBGZ(TZ) rotational speed (RIVU^) of zone IZ about Z-axis 

CBVX(IZ) translational speed of zone IZ in X-axis direction 

CBVY(IZ) translational speed of zone IZ in Y-axis direction 

CBVZ(IZ) translational speed of zone IZ in Z-axis direction 

Record Group #4: Identifies the zonal interface matching indices. (This group input for 

FDNS only, not radiation code.) 

Format: 

NNBC, IZB1, IZF1, IJZ1, UZ2, JKZ1, JKZ2, 

IZB2, IZF2, UZ1, UZ2, JKZ1, JKZ2, 

< (2*IZFACE data lines) 

Definition: 

NNB IZFACE counter (not used in the code) 

IZB1 zonal index of interface plane #1 
IZF1 interface plane identifier for plane #1 
1 : I = I-max (or East) 

2: I = 1 (or West) 

3: J = J-max (or North) 

4: J = 1 (or South) 

5: K = K-max (or Top) 

6: K = 1 (or Bottom) 

IZB2 zonal index of interface plane #2 
IZF2 interface plane identifier for plane #2 

UZ1 the starting point of the first running index on the interface plane 
UZ2 the ending point of the first running index on the interface plane 
JKZ1 the starting point of the second running index on the interface 
plane 

Example: If IZF1 or IZF2 is either 1 or 2 then UZ1 and UZ2 are the indices 

in J-direction and JKZ1 and JKZ2 are the indices in K-direction. 
If IZF1 or IZF2 is either 3 or 4 then UZ1 and UZ2 are the indices 
in I-direction and JKZ1 and JKZ2 are the indices in K-direction. 
If IZF1 or IZF2 is either 5 or 6 then IJZ1 and IJZ2 are 
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the indices in I-direction and JKZ1 and JKZ2 are the indices in J- 
direction. 

Notice: The interface patching surface indices for planes #1 and #2 (i.e. 

UZ1, UZ2, JKZ1, JKZ2) must have consistent running order. 


Record Group #5: 


Specifies flow boundaries (inlet, outlet, symmetry). 
Format: 


IBCZON, IDBC, ITYBC, UBB, UBS, UBT, IKBS, IKBT, 

< (IBND data lines) 


Definition: 

IBCZON 

IDBC 


ITYBC 


UBB 

UBS, UBT 
JKBS,JKBT 


ITYBC 


zonal index for the flow boundary 
boundary facing index 
1: I = I-max (or East) 

2: I = 1 (or West) 

3: J = J-max (or North) 

4: J = 1 (or South) 

5: K = K-max (or Top) 

6: K = 1 (or Bottom) 
identifies boundary type 
-2: inlet fixing everything except pressure 

-1: inlet fixing mass flow rates (e.g. solid fuel blowing 

surfaces) 

0: inlet fixing everything (e.g. supersonic) 

1 : inlet fixing total pressure (compressible flow only) 

2: outlet boundary 

3: symmetry plane (can also be used for slip wall 

boundary conditions) 

I, J or K location (depends on IDBC) of the boundary 
boundary starting and ending indices (for I or J) 
boundary starting and ending indices (for J or K) 

1 : I = I-max (or East) 

2: I = 1 (or West) 

3: J = J-max (or North) 

4: J = 1 (or South) 

5: K = K-max (or Top) 

6: K = 1 (or Bottom) 
identifies boundary type 
-2: inlet fixing everything except pressure 

-1: inlet fixing mass flow rates (e.g. solid fuel blowing 

surfaces) 

0: inlet fixing everything (e.g. supersonic) 

1 : inlet fixing total pressure (compressible flow only) 

2: outlet boundary 


2-43 



SECA-FR-94-18 


3: symmetry plane (can also be used for slip wall 

boundary conditions) 

UBB I, J or K location (depends on IDBC) of the boundary 

IJBS,UBT boundary starting and ending indices (for I or J) 

JKBS.JKBT boundary starting and ending indices (for J or K) 

Record Group #6: Specifies wall block indices. 

Format: 

IWBZON, LI ,L2,M1,M2,N1 ,N2,IWTM,HQDOX,IWALL,DENNX ( VISWX (FDNS) 
IWBZON, LI ,L2,M1 ,M2,N1 ,N2 (TOARAD) 

< (ID data lines) 

Definition: 

IWBZON zonal index for the wall block 

LI, L2 starting and ending indices in the I-direction 

Ml, M2 starting and ending indices in the J-direction 

Nl, N2 starting and ending indices in the K-direction 

IWTM solid-wall thermal boundary condition options 

-1: for fixed wall-temperature 

1: for heat-flux (=HQDOX) b.c. 

HQDOX non-dimensional wall heat flux when IWTM = 1, positive 
from wall to fluid. Normalization for Q : 

SI Units = Q /(p W fU rcf Cp ref T ref ) 

English Units = Q /(32. 174p ren U ren Cp ren T ren ) 

IWALL solid wall heat conduction option 

0: to deactivate; 1: to activate 

DENNX non-dimensional solid wall density 

= (wall-density)/(den-ref) 

VISWX non-dimensional solid wall thermal conductivity 

= k / (x-ref)/(den-ref)/(u-ref)/(Cp-ref) 

Record Group #7: Specifies the singularity lines. (FDNS only) 

Format: 

ISNZON, ISNBC, ISNAX, ISNBS, ISNBT, 

< (ISNGL data lines) 

Definition: 

ISNZON zonal index for the singularity lines 
ISNBC singularity line boundary facing index 

1 : I = I-max (or East) 

2: I = 1 (or West) 

3: J = J-max (or North) 

4: J = 1 (or South) 

5: K = K-max (or Top) 

6: K = 1 (or bottom) 

♦When IWALL = 1 is selected, the program will set IWTM = -1, since this is a correct 
combination. 
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ISNAX orientation of the singularity line axis for example: 
on I-J plane (ISNBC = 5 or 6) 

ISNAX = 1 for I-axis 
ISNAX = 2 for J-axis 
on J-K plane (ISNBC = 1 or 2) 

ISNAX = 1 for J-axis 
ISNAX = 2 for K-axis 
on K-I plane (ISNBC = 3 or 4) 

ISNAX = 1 for I-axis 
ISNAX = 2 for K-axis 

ISNBS, ISNBT starting and ending indices along ISNAX 


Record Group #8: I/O parameters and problem control parameters. (FDNS only) 

Format: 


IDATA, IGEO, ITT, ITPNT, ICOUP, NLIMT, IAX, ICYC, 

< 


Definition: 


(one data line) 


IDATA 


IGEO 


ITT 

ITPNT 

ICOUP 

NLIMT 

IAX 

ICYC 


restart options 

IDATA = 1 for regular restart runs. Restart grid and flow files 
fort. 12 and fort. 13 must be made available. 

IDATA = 2 for example start run. Initial grid and flow data must 
be made available in the fexmpOl include file, 
geometry parameter (for user applications) 

IGEO = 1 is specifically for problems without inlets 
and outlets (e.g. cavity flows) 

IGEO = 9 is reserved for 3-D pump or turbine type 
applications (with ICYC =3) 

IGEO = 19 is reserved for linear cascades applications 
number of time steps limit 

the frequency on printing out solutions (through files fort. 22, 
fort. 23, fort.91, fort.92 and fort.93) 

number of pressure correctors (typically 1 for steady-state 
applications and 3-6 for transient or rough initial start applications) 
typically 1; 

0: for printing out the initial or restart files without going 

through solution procedures 
1: for 2-D planar or 3-D flows 

2: for 2-D axisymmetric flow problems 

cyclic or periodic boundary conditions identifier 
Currently, only ICYC = 3 is active for turbomachinery 
applications. 
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Record Group #9: Time-step size, upwind schemes and time-marching scheme selections 

(FDNS only). 

Format: 

DTT, IREC, REC, THETA, BETAP, IEXX, PRAT, 

< (one data line) 

Definition: 

DTT non-dimensional time step size, DT*Uref/Xref 
IREC selects upwind scheme options 

0 : for second-order upwind scheme 

1 : for third-order upwind scheme 

2 : for second-order central scheme 

REC upwind damping parameter (0. 1 recommended) 

0.0 for second-order accuracy 
1 . 0 : for first-order upwind scheme 

THETA time-marching scheme 0 parameter 
1 . 0 : for steady-state applications 

.99: for implicit-Euler transient applications 

0.5 : for Crank-Nicholson second-order accurate transient 

applications 

BETAP pressure updating under-relaxation parameter typically 1.0; 

small values can be used to reduce the amount on pressure 
corrections for rough start initial runs 
IEXX outlet extrapolation parameter for scalar quantities 
1 : for zero-gradient extrapolation 

2 : for linear extrapolation 

PRAT specifies outlet boundary condition options 
- 1 . 0 : for supersonic outlet b. c. 

0 . 0 : for outlet mass conservation b. c. 

>0.0: for outlet fix pressure b. c. The outlet pressure reference 
point (IPEX, JPEX) is used here. Pressure at this point is 
maintained at a value of PRAT*PPCN. Where PPCN = 
I/ 7 M 2 

Record Group #10: Specifies inlet, outlet pressure points and data monitoring point (FDNS 

only). 

Format: 

IPC, JPC, IPEX, JPEX, IMN, JMN, 

< (one data line) 

Definition: 

IPC, JPC flowfield reference point 

IPC: local grid index in zone JPC (not the global grid 

index) 


2-46 



SECA-FR-94-18 


IPEX, JPEX outlet pressure reference point (same way of indexing as 
IPC, JPC) 

IMN, JMN solution monitoring point (same way of indexing as IPC, 
JPC) 

Record Group #11: Gives reference viscosity, Mach number and options of turbulence models 

(FDNS only). 

Format: 

VISC, IG, ITURB, AMC, GAMA CBE, CBH, EREXT, 

< (one data ii ne ) 

Definition: 

VISC non-dimensional fluid viscosity = l/(Reynolds number) 

= vis-ref/ (den-ref)/u-ref)/(x-ref) 

IG = 1 : for laminar flow option 
= 2 : for turbulent flow option 
ITURB for turbulence model selection 

1: for standard high-Re k-e model 

2: for extended high-Re k-e model 

3: for L-B low-Re k-e model 

4: for H-G low-Re k-e model 

AMC reference Mach number, =(u-ref)/(ref. sound speed) 

GAMA reference specific heat ratio 

CBE non-dimensional buoyancy force parameter = Gr/Re**2, where 
Gr stands for the Grashoff number and Re is the flow Reynolds 
number 

CBH used to activate compressibility corrections for the k-e turbulence 
models 

= - 1 . 0 : for k-corrected model 

= - 2 . 0 : for e-corrected model 

< -3.0: for t-corrected model where Cjfr/T^)' 1 '. 7 = (3- 

CBH) 

EREXT convergence criterion (typically 5.0E-04 for steady-state solutions) 


Record Group #12: Specifies number of zonal iterations in the matrix solver when 

INFACE is used for overlaid grid zonal interface interpolations and 
indicates orthogonal or non-orthogonal grid options (FDNS only). 
Format: 

ISWU, ISWP, ISWK, ISKEW, 

< (one data ji n e) 

Definition: 

ISWU number of iterations for the overlaid zonal boundaries for the 
momentum and energy equations 
ISWU <90: using point implicit matrix solver, 
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ISWP 


ISWK 


ISKEW 


ISWU>90: using conjugate gradient matrix solver with a 

convergence criteria that the residual has to be 
reduced by (ISWU-90) order 

number of iterations for the overlaid zonal boundaries for the 

pressure correction equations 

ISWP <90: using point implicit matrix solver, 

ISWP >90: using conjugate gradient matrix solver with a 

convergence criteria that the residual has to be 
reduced by (ISWU-90) order 

number of iterations for the overlaid zonal boundaries for the 
scalar equations 

ISWK < 90: using point implicit matrix solver, 

ISWK >90: using conjugate gradient matrix solver with a 
convergence criteria that the residual has to be 
reduced by (ISWU-90) order 
non-orthogonal grid viscous flux option 
0: for orthogonal grid 

1: for non-orthogonal grid 


-(one data line) 


Record Group #13: Specifies which equations are to be solved (FDNS only.) 

Format: 

INSO(IEQ): 

U, V, W, TM, DK, DE, 7, 8, 9, VS, FM, SP, 

< 

Definition: (0 to deactivate; 1 to activate) 

U, V, W for the momentum equations 

TM for the energy equation 

DK, DE for the turbulence model 

7, 8, 9 not used 

VS for updating the turbulence eddy viscosity 

FM for the species mass-fraction equations 

SP for calculating the gas thermal properties, and selecting 

various treatment for species production term. 


= 1 

= 11 or 12 

= 21 or 22 

= 31 or 32 

= 33 
> 100 


explicit chemistry model (penalty function) 
implicit chemistry model (1“ or 2 nd -order) 
with psudo-time step size 
implicit chemistry model (1“ or 2 nd -order) 
with real time step size 
implicit chemistry model (1“ or 2“ d -order) 
with time integration to flow time step size 
4th-order PARASOL 

equilibrium plus (SP-100) global finite rate 
chemistry models 
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Record Group #14: Specifies wall radiosity data for IDA method. Input only for IDA 

method. (IDARAD only) Input as many Record 14’s as IDRW. 
IORDR, IDBR, IRADW 

IORDER indicates wall order from which radiosity is to be calculated 
(coincides with ityp order in RAVFAC input file: 
NOZZRAD . INP) 

IDBR = 1 for an IDBC (open) boundary (input, output or 

symmetry); 

= 2 for an ID (solid) boundary 

IRADW boundary index (l:east-I; 2:west-I; 3:north-J; 4:south-J; 
5:top-I; 6:bottom-K) 

Record Group #15: Specifies number of gas species and reactions, and gives the reference 

conditions 

Format: 

NGAS, NREACT, IUNIT, DENREF, UREF, TREF, XREF, (FDNS) 
NGAS,IUNIT,DENREF,UREF,TREF,XREF (IDARAD) 

< ( one data line) 


Definition: 

NGAS 


NREACT 


IUNIT 

DENREF 

UREF 

TREF 

XREF 


number of gas species CEC tables to be read 

= 0: for ideal gas run 

> 0: for CEC real gas run 

=-l: for LOX NBS-table property option 

(Check subroutine INIT for hard-wired LOX 
initial conditions) 
number on reaction steps to be used 
= 0: for non-reacting flow 

>0: for finite-rate reacting flow 

= 1: for Si-unit reference conditions 

= 2: for English-unit reference conditions 

reference density (in kg/m 3 or slug/ft 3 ) 
reference velocity (m/sec or ft/sec) 
reference temperature (°K or °R) 
reference length (m or ft) 


Record Group #16: Include the CEC thermodynamics data here 

Format: 

Name, Molecular Weight, Coefficients (7 x 2) 

< (4*NGAS lines) 

FDNS reads in the data in CEC format. 
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Record Group #17: Specifies the finite-rate reacting steps (FDNS only) 

Format: 

REACTING: Species Names, N = 1, NGAS (this is a title line) 
IREACT, A, B, E/RT, ITHIRD, IGLOB 
(STOCEF(N, IREACT), N=1,NGAS) 

(STOCEG(N, IREACT), N=1,NGAS) — If IGLOB = 2 

< (NREACT sets) 

Definition: 

IREACT reaction step counter 

A reaction rate leading constant 

B reaction rate temperature exponent 

E/RT reaction rate activation energy constant 

ITHIRD third-body reaction indicator 
0: deactivated 

N: for using the N-th species as third body 

999: for global (every species) third-body 


Record Group #18: provides particle input control 
Format: 

IDPTCL, IPREAD 

< (1 data line) 

Definition: 


IDPTCL number on particle sizes initial condition input lines 
0: to deactivate particulate phase option 

IPREAD 1 for reading in particle data (fort. 14) from upstream domain (this allows 
transferring the outlet particle data from the upstream domaine 
solutions to the inlet boundary for succeeding domain computations 
—especially useful for multi-phase rocket plume simulations) 
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Record Group #19: for reading in particle initial conditions (for steady-state runs only) (FDNS 

only) 

Format: 

IPTZ0N,IDBCPT,LPTCL1,LPTCL2,MPTCL1,MPTCL2,NPTCL1,NPTCL2, 
ITPTCL,DDPTCL,DNPTCL, WDMASS ,UUPTCL,HTPTCL 

< (2*IDPTCL data lines) 

Definition: 

IPTZON zonal index for the particle initial position 
IDBCPT I-, J- or K-plane identifier 

1 : for I-plane (plane normal to I lines) 

2: for J-plane (plane normal to J lines) 

3: for K-plane (plane normal to K lines) 

LPTCL1,LPTCL2 I-interval for the particle initial position 

MPTCL1,MPTCL2 J-interval for the particle initial position 

NPTCL1,NPTCL2 K-interval for the particle initial position 

ITPTCL number of particle groups (trajectories) starting from each grid cell 

DDPTCL particle diameter in /xm 

DNPTCL particle density in lbm/ft**3 

WDMASS particle mass flow rate for the current particle group 

UUPTCL particle/gas velocity ratio at the initial positions 

HTPTCL particle initial enthalpy in ft**2/sec**2 
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(6) Execute the RADO program. This code generates the fort. 15 and fort. 16 initial 
guesses for radiant intensities for the gas and particles at each flowfield point for 
each of the gas and particle bands that are specified by the user in the 
BLKDAT_*.f routine. 

(7) Execute the ID ARAD program. Output of the code consists of a fort. 6 file that 
tracks the progress of the code through each of the bands (gas and particle) as the 
program executes and can provide the user information necessary to correct any 
errors that are encountered during execution. The actual radiation output in the 
form radiation heat fluxes at each wall (or boundary point) are found in fort. 37 
and fort. 67. Fort. 37 gives the location of the wall point along with the difference 
in incident radiation to the point and radiation intensity at the point due to the 
medium at the point. Fort. 67 provides the axial location of the point and the 
incident radiation at the point (wall or boundary). For ODA cases, the radiation 
intensities (BTU/ft 2 /sec) are output for all wall points. In the case of an IDA 
calculation, only those points that are specified as radiating boundary points have 
non -zero heat fluxes. 

Three data files that must reside in the directory in which the IDARAD and RADO codes 
are being executed in are: SIGACL.DAT, SIGSCL.DAT and koch2.pm. SIGACL.DAT 
contains the absorption coefficient data for A1 2 0 3 . SIGSCL.DAT contains the scattering 
coefficient data for A1 2 0 3 and koch2.pm contains the extinction coefficients for A1 2 0 3 . 
Particulates other than A1 2 0 3 would require regeneration of the files to describe the properties 
of the particular particulate. 

Additional notes on the options available to the IDARAD code are shown in Table 2. 12. 
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Table 2.7 Listing of fort. 11 ID ARAD Sample Input File 


IDIM, 

r> 

(2-D asrm48-5 nozzle. 

PARTICULATE TWO- 

PHASE 

FLOW) 

IZON , 

IBND, 

ID, 

I RAD, 

IDRW , 

IFL13 




1, 


1 , 

4, 

4, 

1 




IZT, 

JZT , 

KZT , 







179, 

81, 

1. 







IBCZON . 

IDBC, 

ITYBC , 

IJBB. 

IJBS. 

I JBT , I KBS, 

I 

KBT , 

1, 

1, 

o 

<• 

179. 

1. 

81, 

1, 


1. 

1, 

fl 

1, 

1, 

1, 

81, 

1, 


1, 

1, 

4, 

3 , 

1, 

1, 

179, 

1, 


1, 

IWBZOM , 

LI, 

L2, 

Ml, 

M2, 

HI, 

N2 , 



1, 

1. 

179, 

81, 

81, 

1, 

1, 



I ORDER, 

IDBR , 

IRADW, 







1, 

1, 


! rad 

segmen t 

1 , boundary 

1, 

west i 

**> 
*- , 

1. 

, 

! rad 

segmen t 

2, boundary 


west 1 

3 , 

n 

, 

3 , 

f rad 

segmen t 

1, boundary 

*7* 

north 

4, 


3 , 

! rad 

segmen t 

2, boundary 

n 

n 

north 

NGAS, 

IUWIT, 

DREE(SLG), UREF(FZS), 

TREF(R) , 

XREF 

(FT) , 

12, 

2, 

5 .299 IE— 

03, 116-68630, 

6319.80, 

1 .00000000 , 


H20 

0.26340654E+01 
-0 . 29876258E+05 
—0 . 48670872E-08 
02 

0.36122139E+01 
-0.1197815IE+04 
—0 .9818910 IE-08 
H2 

0.30558124E+01 


300.000 5000.000 18.01520 

0 . 31 121899E-02-0 . 90278451E-06 0 . 12673054E-09-0 .69164734E-14 
0.70823874E+01 0 . 41675563E+01-0. 18106868E-02 0 . 59450877E-05 
0 . 1 5284144E-1 1-0 . 30289547E+05-0 . 73087996E+00 

300.000 5000.000 31-99880 

0 . 74853166E— 03— 0 . 19320646E-06 0 . 33749007E-10-0 . 23907374E-14 
0 . 36703308E+01 0 . 37837136E+01-0 . 30233634E-02 0 . 99492754E-05 
0 . 33031826E— 1 1— 0 . 10638107E+04 0.36416345E+01 

300-000 5000.000 2.01580 

0 . 59740403E— 03— 0 . 16747471E-08-0 . 21247544E-10 0 . 25195486E-14 
— 0 . 86168475E+03— 0 . 17207073E+01 0 . 2943232SE+01 0 .34815508E-02-0. 77713821E-05 
0. 74997493E-08-0.25203379E-1 1-0. 976954 10E+03-0. 13186136E+01 
0 300.000 5000.000 15.99940 

0.25342960E+01— 0. 12478170E— 04-0 . 12562724E-07 0 . 69029860E-1 1-0 . 63797098E-1 5 
0.29231107E+05 0 . 4962B592E+01 0 . 3030940 1E+0 1-0 . 22525853E-02 0 . 39824 540E-05 
— 0 . 3260492 1 E-08 0 . 10 1 52035E-1 1 0 . 29136525E+05 0 . 26099341E+01 

H 300.000 5000.000 1.00790 

0.25000000E+01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 
0 . 2 547439 1E+0 5—0 - 45989841E+00 0 . 25000000E+01 0 . 00000000E+00 0 . 00000000E+00 
0.00000000E +-00 0 - 00000000E+00 0 . 2 547439 IE +0 5—0 . 4 598984 1E+00 

300.000 5000.000 17.00730 

0 . 10005879E— 02— 0 . 22048808E— 06 0 . 20191288E-10-0 . 39409830E-1 5 
0. 55566425E+01 0 . 38737299E+01-0 . 13393773E-02 0 . 1634S351E-05 
0 . 4182697 5E— 13 0 . 35802349E+04 0 . 34202406E+00 

300.000 5000.000 28.01040 

0 . 14891390E— 02-0 . 57899683 E-06 0 . 10364577E-09-0 . 69353550E-14 
0 . 634791 56E+01 0 . 371 00928E+01-0 . 16190964E-02 0 . 36923593E-05 
0.23953344E-12— 0. 14356310E+05 0 . 29555352E+01 

300.000 5000.000 44.00980 

0 . 3098 17 19E— 02-0 - 1239257 IE— 05 0 . 22741325E-09-0 . 1 5525955E-13 
-0 . 48961 44 1E+0 5—0 . 98635983E+00 0 . 24007797E+01 0 . 87350961E-02-0 . 66070879E-05 
0.20021362E-08 0 . 63274039E-1 5-0 . 48377527E+05 0 . 96951456E+01 
CL?7? 300.000 5000.000 35.45300 

0 . 29537796E+01— 0 . 40792712E— 03 0 . 1 5288342E-06-0 .26384345E-10 0 . 17206581E-14 
0.13695677E+05 0 . 30667325E+01 0 . 2077428 1E+01 0 . 29487169E-02-0 . 43919732E-05 
0 . 24499776E-08— 0 . 41007685E— 12 0 . 13871 928E+05 0 . 73136343E+01 
CL277 300.000 5000.000 70.90600 

0 . 430778 14E+01 0 .31 182816E-03-0 . 16310807E-06 0 . 4451 1913E-10-0 . 43057753E-14 
— 0 . 13458251E+04 0 . 20666684E+01 0 . 31316886E+01 0 . 48997877E-02-0 . 6941 1463E-05 
0 . 4478564 1 E-08— 0 . 10621859E-1 1-0 . 10979696E+04 0 . 77833424E+01 
HCL77 300.000 5000.000 36.46100 

0 . 27665884E+01 0 . 14381383E-02-0 . 46993000E-06 0 . 73499408E-10-0 . 43731 106E-14 
— 0. 11917468E+05 0 . 64583540E+01 0. 35248171E+01 0. 29984862E-04-0. 86221891E-06 
0 . 2097972 1 E-08— 0 . 98658191E-12-0 . 121 50509E+05 0 . 23957713E+01 
N2 300.000 5000.000 28.01340 

0 . 28 532898E+0 1 0 . 1 6022 128E-02-0 . 6293689 IE-06 0.11441 022E-09-0 . 780 5 7466E- 1 4 
— 0 . 89008093E+03 0 . 63964896E+01 0. 37044 177E+0 1-0 . 14218753E-02 0 . 28670393E-05 
-0 . 1 202888 5 E-08— 0 . 13954677E-13-0 . 1064079 5E+04 0 . 22336285E+01 


OH 

0 .2889781 5E+01 
0 . 388 5704 1E+04 
-0.52133636E-09 
CO 

0.29840696E+01 
—0. 14245228E+05 
— 0.203 1967 5E-08 
C02 

0 . 44608040E+01 
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Table 2.8 Description of IDA Input File NOZZRAD.INP* 


Record 1: Header Format:: 80A1 

Variable Column 

HEAD 1-80 Problem description 


Record 2: 


Variable 

nstart 1 Restart control flag 

nt 0 Output tape for restart 

nvfcal 1 View factor calculation technique option (also used to calculate 

solid angles). The solid angle computations are coded for the 
finite difference technique only. Therefore, nvfcal is set to 1 in the 
ravfac_sa subroutine. 

norm 0 View factor normalization option 

rmax 0.0 Maximum area- to-di stance ratio 

n Prt 0 Immediate output control 

nfe 0 Element override option 

nfs 0 Shading override option 

ntvf 1 View factor output tape 

Record 2: Format: Free 

IPLANE 

Variable Value Description 

IPLANE N Number of 3D planes (suggest 2) 


Record 4: Format: Free 

IDRW 

Variable Value Description 

IDRW N Number of radiating boundary segments 

♦Each input record has a header card associated with it. See Table 2.9 for sample case. 


Format: Free 

nstart, nt, nvfcal, norm, rmax, nprt, nfe, nfs, ntvf 
Value Description 
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Table 2.8 Description of IDA Input File NOZZRAD.INP (Continued) 


Record 5: Format: Free 

IT YPE(I) , IT YPE(2) ITYPE(IDRW) 

Variable Value Description 

(ITYPE(I),I=1,IDRW) 2 Circular plane 

5 Cone 

This record specifies the geometry used to describe each IDRW segment. A negative 
sign ahead of the variable is used to specify an outside surface. 


Record(s) 6: Format: Free 

ILV1, ILV2, LQ, MV1, MV2, MV3 

Record 6 is input for each IDRW surface which describes the flowfield grid indices that 
describe the inlet and wall joints that are to be treated as radiating surfaces. 


Variable 

Value 

Description 

ILV1 

N 

Beginning i index 

ILV2 

N 

Ending i index 

LQ 

N 

i increment 

MV1 

N 

Beginning j index 

MV2 

N 

Ending j index 

MV3 

N 

j increment 


Record 7: Format: Free 

ILOOP, INDIV, INRB 


This record inputs the 
of three input variables. 

Variable Value 

ILOOP N 

INDIV N 

INRB N 


number of radiation source flowfield segments or points for each 


Description 

number of do loops used to input flowfield points 
number of individual points to input as flowfield points 
number of non-radiating boundary flowfield do loop’s 
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Table 2.8 Description of IDA Input File NOZZRAD.INP (Continued) 


Record(s) 8: Format: Free 

INI, IN2, INQ, JN1, JN2, JNQ 

Input one record 8 for each ILOOP zones 


Variable 

Value 

Description 

INI 

N 

Beginning i index for iloop flowfield zones 

IN2 

N 

Ending i index for iloop flowfield zones 

INQ 

N 

i increment 

JN1 

N 

Beginning j index for iloop flowfield zones 

JN2 

N 

Ending j index for iloop flowfield zones 

JNQ 

N 

j increment 

Record(s) 9: 

Format: 

Free 

IXNODE, IYNODE 

Input one record 9 for each specified INDIV flowfield point 

Variable 

Value 

Description 

IXNODE 

N 

i point index 

IYNODE 

N 

j point index 

Record(s) 10: 

Format: 

Free 

INI, IN2, JN1, JN2, JNQ 

Input one record 10 for each non-radiating boundary point (INRB) loop 

Variable 

Value 

Description 

INI 

N 

Beginning i index for points 

IN2 

N 

Ending i index for points 

INQ 

N 

i increment 

JN1 

N 

Beginning j index for points 

JN2 

N 

Ending j index for points 

JNQ 

N 

j increment 
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Table 2.8 Description of IDA Input File NOZZRAD.INP (Continued) 

General notes on inputting NOZZRAD.INP 

The input files contain headers for each input record. The headers must be input as 
shown in the sample input case (Table 2.9). 

Two ranges of interpolation domains are required in the input (inlets/walls and flowfield.) 

User must make sure that the comer and edge points of each interpolation domain are specified 
as either radiating boundary points or flowfield points. This ensures that interpolation of 
unknown incident radiation points in INTRP.f will be bounded by known values. These points 
are indicated in the input file ’NOZZRAD.INP’. 

a. The radiating boundary points must be bounded by known values. This is accounted for 
by providing the node values associated with (Ivl,lv2,lq;mvl,mv2,mq), and must 
include the endpoints and comer points of the radiating boundaries. 

b. The flowfield points (not including radiating boundary points) must be bounded by known 
values. This is accounted for by providing the node values associated with iloop, indiv 
and inrb. The limits of the incident radiation nodes should extend one point off the 
radiation boundary nodes. 

Incident RAD or flow nodes: 

Do not input overlapping or duplicate points; i.e., each point specified through the iloop, 
indiv and inrb (non-radiating boundary) parameters must be unique (see subroutine 
VFINIT.f for more description). 

Radiating boundary nodes: 

Boundary points specified through (Ivl,lv2,lq;mvl,mv2,mq) can overlap, especially 
when merging two segments together. 

Can input segments in any order. Must specify at least two points per radiating boundary 
segment (can specify open flow (IDBC) boundaries as radiating boundaries). 

Do not overlap radiating boundary points with incident flow points - this will result in incorrect 
calculation of interpolation indices in SORT.f. 

It is important that the incident radiation source flow nodes input in NOZZRAD.INP bound all 
of the interior nodes, not including the radiating boundaries (no error message is generated if 
this procedure is not followed). 
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Table 2.8 Description of IDA Input File NOZZRAD.INP (Continued) 


If iloop range, including step values, does not extend to all non-radiating boundaries, must 
include points via indiv or inrb input regions. Otherwise, SORT.f will not run correctly. 

When running the IDA case within a portion of the total grid domain, the rule stated above must 
be followed (for 2D); user must supply radiating boundary such that all comer and edge points 
are specified, and incident radiation points such that all comer and edge points are specified; 
all of the comer and edge points together must form a square region. 

Interpolation indices are determined only for J-lines that have more than one known point. 

To get the best interpolation between known IDA points, and to avoid interpolation within 
skewed cells, place the radiating boundary points in the same J-line as the flowfield node points. 
Otherwise, there may result only one point in a given J-line (that corresponding to the radiating 
wall point). 

When splitting a radiating boundary into more than one segment (Ivl,lv2,lq;mvl,mv2,mq), user 
needs to group together all segments associated with the same wall in the RAD wall input 
portion of the NOZZRAD.INP file. 

The differential increments along a line-of-sight in the Z (DSZ) and R (DSR) directions are 
currently set at 0.1 ft in subroutine OPDIST.f. 


Zones 

The radiation coding is written for a single grid zone only (IZON=l). However, there 
can be more than one radiation zone (IRADZN) through the inputs corresponding to 
iloop, indiv and inrb in file NOZZRAD.INP. 

Wall Boundary 

The radiation portion of code assumes following wall boundary location: A single nozzle 
wall located on north boundary. 


VFINIT.f notes 

Initializes the RAVFAC input file for a nozzle 

Can choose a number of radiating wall points through the increment in the axial (LQ) and radial 
(MQ) directions (required input for view factor and solid angle runs) 

Can also choose a number of flowfield grid nodes (required input for solid angle run) 
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Table 2.8 Description of IDA Input File NOZZRAD.INP (Continued) 

The order of accuracy in calculating the view factors and the solid angles at each point can also 
be adjusted through the factor ibe(# of elements in beta dir) and 

ige(# of elements in gamma dir). 

vfinit reads in NOZZRAD.INP and generates the file: 

A. IN.DAT_VF when iflow = 0 

View Factor Pre-processor - 

view factor computations from wall point to wall point: 
used in determining radiosity RADOSW in IDA method. 

B. IN.DAT_SA when iflow > 0 

Solid Angle Pre-processor - 

solid angle computations from wall point to flowfield node 
used in determining incident radiation RIO in IDA method. 

iloop - number of do-loop input ’zones’ 

indiv - number of individual point input ’zones’ 

inrb - number of non-radiating boundaries 

1 . iflow is set to 1 for wall points to flowfield points by inputting through a do-loop 
(iloop >0). 

The order of input is: do loop for radial dir. 

do loop for axial dir. 

This will result in 2 grid regions if all flowfield nodes are chosen: 

FDNS flowfield grid, IDA grid = SIRRM grid, 
or in 3 grid regions if some flowfield nodes are chosen: 

FDNS flowfield grid, SIRRM grid and IDA grid. 

To get incident radiation at all points in flowfield with IDA method, program user must 
supply some nodes on all boundaries in the input (through iloop, indiv and/or inrb). This 
will provide values on all boundaries so that interpolation can be effective. 

If a do-loop ("‘iloop > 0) with an increment in the i-node does not allow the max i-node 
to end on a boundary, then must input some boundary points individually or must specify 
points on a non-radiating boundary. 

(If iflow "ge" 1, points chosen lie in axisymmetric plane) 

iplane - number of planes within total angle (g2-gl) 

The circumferential extent of the wall surfaces is taken into account by the angle 
parameters gl and g2. 
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* 


** 


Table 2.9 Listing of NOZZRAD.INP Sample Case 


NOZZLE VIEW FACTORS 

nstart nt nvfacl norm rmax nprt n-fe nsf ntvf (nvfacl must = 1 ) 
1 0 1 0 0 . 0001 
H o-f planes in 360 degrees (i plane) 


— RADIATING WALL OR RADIOSITY POINTS — 
number of radiating walls (idrw) 

4 


type o-f wall (ityp) C 2 =circular plate, 5 =cone; (-) 

2 2-5-5 

indices for radiating wall 1 ( lvl , 1 v 2 , lq ;mvl ,mv 2 ,mq ) 
111 1 77 19 

indices for radiating wall 2 ( lvl , lv 2 , lq ;mvl ,mv 2 ,mq ) 
1 1 1 77 81 A 

indices for radiating wall 3 ( 1 vl , 1 v2 , lq ;mvl ,mv2.mq ) 
1 177 44 81 81 1 

indices for radiating wall 4 ( lvl , lv2, lq ;mvl ,mv2.mq ) 
177 179 2 81 81 1 


for outside 

( 5 )* 

( 2 ) 

(5) 

( 2 ) 


— RADIATION SOURCE FLOWFIELD OR INCIDENT RADIATION POINTS — 
# of do— loop zones (Hoop)* individual points (indiv)* 
and non— rad boundaries (inrb) 

o'?'' 

4 - 4 . iL 


surface] 


a. clo-loop poin ts: in 1 , in 2 „ inq ; j n 1 j, j n 2 , ;i nq 

45 177 44 1 77 19 

45 177 44 80 80 1 

b. ixnode* iynode indices for flowfield points: 

2 80 

179 80 

c. non-radiating boundary points: in 1 „ in 2 „ inq 5 ini * in 2 „ inq 

2 2 1 1 77 19 

179 179 1 1 77 19 


(20) ** 

(4) 

( 1 ) 

( 1 ) 

(5) 

(5) 


The variable IWR (# of radiating boundary points) in BLKDAT_*.f must be set based 
on the wall points selected in this input file which is 14 for this case. 

The variable IWN (# of radiating flowfield points) in BLKDAT *.f must be set based 

on the flowfield points selected in this input file which is 36 for this case. 
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PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 


Table 2.10 Description of radOl Include File 

(IIQMAX = 14500, IWP = 201) 

!max grid points and max wall points in 2D 
(NSPM = 12, ISPMAX = IIQMAX) 

!max species and switch that lets code know that every point in domain 
requires species calculations (set ispmax to 1 otherwise) 

(NPMAX = 1, UKPMX = IIQMAX) 

! number types of particles and switch that lets code know that every point 
in domain requires particle calculations (set ijkpmx to 1 otherwise) 
(IMAP = 60) 

!60 SIRRM map input nodes (for particle radiation) 

(IPLMAX = 2, IDRWMX = 4, IDRNMX = 6) 

! Two 3D revolved planes for IDA radiation 
(IWR = 14, IWN = 36) 

! number of radiating boundary points and incident radiation nodal points 
in interior for IDA) 

(NGBAND = 4, NPBAND = 5) 

!max number of gas and particle bands or windows (choice of npband is 
tricky when irad=4 since the number of particle band segments the code 
chooses may not be evident. Choose a reasonable value for this variable). 
(IWRP = IWR*IPLMAX, NTBAND = NGBAND + NPBAND) 

Notes on setting variables in RAD01 


Change IIQMAX and IWP in RAD01 for appropriate dimensions, (for multi-species problems 
set NSWPM = number of species and ISPMAX = IIQMAX. Note: Do not change other 
parameters.) 


Set IMAP to the larger of the following values: 

IXSTA (number of x-stations in SIRRM map file) 

IYSTA (number of y-stations in SIRRM map file) 

Set IPLMAX to the number of 3D planes required for radiosity calculations. 

Set IDRWMX to the maximum number of radiating walls in the axisymmetric plane, and 
IDRNMX to the maximum number of radiating source zones in the axisymmetric plane (see 
NOZZRAD . INP file). 

Set IWR to the maximum number of radiating boundary points on any boundary for IDA in the 
axisymmetric plane (if two boundary comers have overlapping point, must count as 2 points). 
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Table 2.10 Description of radOl Include File (Continued) 

Set IWN to the maximum number of flowfield node points for IDA incident radiation in the 
axisymmetric plane. (IWR and IWN must be equal to the values in NOZZRAD.INP that 
correspond to the number of radiating wall or radiosity points (lvl, lv2, lq:mvl,mv2,mq) and 
to the number of radiating source or incident radiation points (iloop;indiv;inrb), respectively). 

Set NGBAND to the dimension of the total number of gas bands, summed over all gas species. 

Set NPBAND to the dimension of the total number of particle bands (or particle window bands 
for I RAD =4), summed over all particle species. For IRAD=4, make NPBAND a little larger. 
This is required since the number of particle bands may be split into additional band segments 
due to gas/particle overlap. 

(ex.l: 

for IRAD=3, NGBAND = 1 1 for all gas bands of H 2 0 and C0 2 ; 

NPBAND = 12 for all particle bands from 0.5 to 6.0 microns.) 

(ex2: 

for IRAD=4 NGBAND = 1 1 for all gas bands of H 2 0 and CO z ; 

NPBAND = 18 for all particle bands from 0.5 to 6.0 microns, 
where the 18 bands are associated with the particle 
window regions (those regions that are not 
associated with the gas/particle overlap regions.)) 
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Table 2.11 Notes in Setting BLKDAT_*.f Variables 

There are several options to choose to run a radiation case, as described in the following fortran files: 

BLKDATO.f (radiation initializing), 

BLKDATU.f (uncoupled radiation code) 

NOTE: The BLKDAT _*.x files differ in the amount of data required to run the specific code in their respective 
subdirectories. Make sure that any data changes made to BLKDAT_0.f is likewise made in BLKDAT U.f. 
Also, recompile the code when BLKDAT *. f has been modified. 

However, most of the values in the block data files above can be left as the default values. 


The most common values which must be changed are the following: 


1 . IAP: 

= 1 : OD A (optically thick region only) 

= 2: Defunct 

= 3: IDA (all optical thicknesses) 


2. For gas cases: 

NSPMS, NSPME - The starting and ending indices for the gas species 
NLAMGS,NLAMGE - The corresponding wavelength indices as listed in The ALAMG data 
statement 


3. For particle cases: 

NPMAXS,NPMAXE - The starting and ending indices for the particle species 
NLAMPS,NLAMPE - The corresponding wavelength indices as listed in The ALAMP data 
statement. 


4. IDSPG 


Indices corresponding to the gas species number as input in the CEC section of fort. 1 1 (those 
corresponding to H 2 0 and or COJ 


Example: To run an IDA case with gas/particle overlap 
NSPM = 1 (starts with H 2 0) 

NSPME = 2 (ends with COj) 

NLAMGS(l) = 3 (starts with 2.7 microns for H 2 0) 

NLAMGS(l) = 4 (ends with 1.87 microns for H 2 0) 

NLAMGS(2) = 4 (starts with 4.3 microns for COj) 

NLAMGS(2) = 5 (ends with 2. 7 microns for COj) 

NPMAXS — 1 (starts with A1 2 0 3 ) 

NPMAXE = 1 (ends with A1 2 0 3 ) 

NLAMPS(l) = 7 (starts with 1.0 micron for A1 2 0 3 ) 

NLAMPS(l) = 9 (ends with 3.0 microns for A1 2 0 3 ) 

SET IVIEW — 0 for ODA cases and 
IVIEW = 1 for IDA cases 

Additional notes on the variables set in BLKDAT*. f can be found in comments contained in BLKDAT*. f 
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Table 2.12 Additional Notes on the Options Available in ID ARAD 


Number of gas-only bands: NOLAMG 

Number of particle-only bands: NOLAMP 


Wavelength sequence (ILT) for ISPECL > 0: 

IRAD = 1 (gas-only): 

ILT sequence corresponds to descending wavelength order 


ILT = 0 

DO ISP=NSPMS,NSPME 
DO IWV=NLAMGS(ISP),NLAMGE(ISP) 

ILT = ILT+1 
END DO 

IRAD = 2 (particle-only): 

ILT sequence corresponds to ascending wavelength order 
ILT = 0 

DO IPA=NPMAXS,NPMAXE 
DO IWV +NLAMPS(IPA),NLAMPE(IPA) 

ILT = ILT+1 
END DO 

IRAD = 3 (non overlapping gas and particles): 

Determines ILT according to IRAD = 1 first, then continues ILT according to IRAD = 2. 


(outer loop) 
(inner loop) 


(outer loop) 
(inner loop) 


IRAD=4 


* gas-only bands 

* Overlapping gas and particles within the total particle band width DLAMP (continuum 
width DLAMP chosen by the user), 

* The remaining particle bands that fall within DLAMP 

* The DLAMP bands that have no overlapping gas bands (referred to as windows): 


a: Determines gas-only and gas/particle overlap sequencing first by ranging through the loops as in 

IRAD = 1 . The wavelengths are in descending order, consistent with the input order for ALAMG 
above. (This descending order must be adhered to as the gas band widths in WIDEMB.f are 
calculated accordingly). 

b. Continues ILT sequencing by finding the remaining particle bands that do not overlap with any 
gas band (window region). These particle wavelengths are in ascending order, with the particle 
band widths calculated as described in i and ii below. (Gas wavelength parameters are computed 
in subroutine RAD0(1) so that they are in ascending wavelength order, making it more convenient 
to check for gas/particle overlapping with the ascending wavelength of the particles, and separating 
out the remaining particle bands from the gas bands). 

i. If there are overlapping gas bands, the total particle band width (DLAMP) is split into 

a number of smaller particle bands (reduced particle band width) which lie between the 
two DLAMP endpoints ALPEND. The particle bands extend from ALPEND to the 
nearest gas band endpoint, and/or between gas band endpoints. 
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Table 2.12 Additional Notes on the Options Available in ID ARAD (Continued) 

ii. If there are no overlapping gas bands within the particle width DLAMP, the particle band 

corresponding to ILT is ALAMP. 

NOTE: The number of particle-only band ’centers’ (ILT sequencing) may not equal the number of particle-only 
band widths. For example, when a gas band overlaps or extends past the ALPEND boundary (gas band 
width varies as the temp., pressure and species mass fraction varies), the reduced particle band width 
would equal zero since there would be no particle gap (See coding in PBAND.f). 

NOTE: For IRAD 4 

a. The portion of the particle band that overlaps with a gas band is accounted for. The effect of overlapping 
particle bands are not subtracted. 

b. If two neighboring particle band widths overlap, there will be a discrepancy in the value of NOLAMP 
(should be decreased by 1). For different press and temp states, NOLAMP could feasibly fluctuate. The 
code requires a constant value for NOLAMP. Until a method is devised that computes a single PI equation 
for all of the particle-only bands together (which would automatically take the varying band width total into 
account), NOLAMP must remain constant. Therefore, if two particle bands overlap, the value for 
DLAMG is set to zero. 

c. Band widths which extend beyond the input values for ALAMP are necessarily cut off at the band end 
points (ALPEND) of the band centers (ALAMP). This shortcut reduces the coding complexity. The result 
is an increase in the total particle-only bands, which may be somewhat erroneous. However, the gas band 
widths are accounted for without shortcuts, since the actual value for the gas band widths are used. If this 
causes concern, the remedy is to choose new (initial) values for ALAMP (and corresponding values for 
ALPEND) which bypass this problem. For example, if the 1.38 band of H 2 0 has a width that extends 
beyond the 1.5 particle band end point, the value 1.5 could be changed to 1.7 to allow the 1.38 band width 
to fall within the particle band search region; or, the value of ALAMP could be kept as is, but ALPEND 
could be adjusted. 

d. For more than one type of solid species (eg, A1 2 0 3 , C(S), etc), the loop (IPA in subroutine RADINO) 
would have to be extended to the other species, in addition to a single species. 

eg.: DO IS = 1,NSPEC 

DO N=l, IRP 

This would be required when choosing the option IRAD=4, where the extinction coefficient (averaged over 
all particle sizes) at a gas band center would be a summation of the SIGEXT of the gas band plus the 
overlapping particle bands. No provision is made in the code for this at present; only a single solid 
species type (A1 2 0 3 ) is allowed (SIGEXT = SIGABS gas + SIGEXT_A1 2 0 3 ) 

Wall Emissivity 

The code is presently set up to allow the same wall emissivity value at evety wall point. Also, the black 
wall option, IBLAKW, must be set to 0 for IDA cases, since the black wall boundary condition case is not 
coded. 

ODA Case 

If a larger number of ODA iterations is required for additional convergence, along with a tighter tolerance 
on the two convergence criteria, these values can be changed in subroutine ODA.f: 

IODAIT = # of iterations (default = 1500) 

TOL1 = Tolerance for average residual (default = l.E-9) 

TOL2 = Tolerance for maximum residual (default = l.E-6) 

Additional description can be found in read.meradul, read.me _rad_u 2 (uncoupled radiation code) 
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2.2.3.2 Installation of IDARAD 

ID ARAD was developed and checked out on a UNIX based IBM RISC work station. 
There is no system specific coding so that the program can be readily compiled on any system 
that has a fortran compiler and sufficient storage. The core storage requirements are dependent 
on the problem being run (i.e., dimensions set in radOl). The radiation initiation code (RADO) 
and the actual IDARAD code (RAD) should be loaded in separate directories. The RADO 
initialization code is contained on the MS-DOS disk RADI. The IDARAD code is contained 
on MS-DOS disk RAD2. Table 2.13 lists the make file that contains the compilation and links 
instructions for the flow initialization code (RADO), using an IBM AIX XL/6000 fortran 
compiler. Table 2.14 contains a listing of the make file for the IDARAD code. Table 2.15 
contains a listing of the functional subroutines and include files that make up the IDARAD code. 

The most efficient way to use the initialization and IDARAD codes is to execute them 
in a separate directory for each problem. In addition to fort.ll, fort.12, fort.13, fort.61 and 
NOZZRAD.INP files that are set up by the user, the optical properties files; koch2.pm, 
SIGACL.DAT and SIGSCL.DAT must also be contained in the working directory. The MS- 
DOS disc RAD3 contains sample case input data files and the optical properties files for running 
the cases whose results are presented in the next section. The fort. 1 1 files are fort. 1 l_ODA and 
fort.ll_IDA for the ODA and IDA cases. 

2.2.4 IDA and ODA Results 

The experiment selected to check out the IDARAD code was the MNASA 48 inch 
contoured ASRM nozzle plume radiation test. A Cycle 2.0 SIRRM map was converted into 
fort.12, fort.13 and fort.61 grid, flowfield and SIRRM mapped files for a 20 x 21 grid. The 
fort.ll, fort.12, fort.13, fort.61 and NOZZRAD.INP files are contained on the MS-DOS disk 
RAD2. ODA and IDA results consist of emissive power at the boundary as a function of axial 
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Table 2.15 Listing of Fortran and Includes Files for the IDARAD Program 


FORTRAN FILES: 

* rinit.f » initialization files 
VFINIT.f RAVFAC.f PATHD.f PATHL.f 
rl-f r2.f radpl.f 

* r prop- f = particle input files 
SIGACL-f SIGSCL.f REFRIN-f 

* rgas.f = gas property files 

ECOEFG.f WIDEBM.f SUHMAl.f SUrtMA2.f ASTAR.f 

* rprtcl.f = particle property files 

ECOEFP-f SIGAV.f REFRAC.f PBAND.f OVRLAP.f 

* rad in t - f = radiation driver files 

RADIN0 . f MONCHR-f EXTINC-f PLANKF.f HTRAD-f 

* roda-f = oda files 

□DA-f TRDIAG.f ODABC-f DLAMB - f PLANKB.f 

* rida-f = ida files 

IDA.f GWALL.f GMEDUM-f OPDIST.f XDABC-f 

INRAD. f 

* rlos.f = line-of-sight files 

LSIGHT.f TAUDR.f SRCIDA.f 

* rwall-f = radiosity files 

RDOSTY-f JWALL.f INVERT. f 

* rintrp.f = interpolation files 

LOGINT.f LININT.f LININR.f CELLIJ.f SORT.f 
INTRP0 . f INTRF'l - f INTRP2.f INTRP3.f 

* BLKDATJJ.f « block, data file 


INCLUDE FILES: 

common block files - 

rad01 ? rad02 , . . - 9 rad24 , 
and alumox -inc, radl.inc 
fortran code include files — 

rad2.inc - error statement check 

init-inc ~ sets additional radiation parameters 
dataiol.inc - reads in radiation files 
fort. IS and fort. 16 
dataio2.inc - outputs radiation files 
fort. 25 and fort. 26 
widebml.inc & widebm2.inc - 

additional fortran code for 
wide band model (WIDEBM.f) 
h2o . in c - fortran code for determining 
wide band parameters for h2o 
co2. inc - fortran code for determining 
wide band parameters for co2 
print_u- inc - heating rate driver 
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MNASA Predictions 



Axial Distance Fron Nozzle Exit (ft) 

Comparison of ODA and IDA Prediction of Plume Emissive Power Distributions 
for the MNASA48 inch Contoured ASRM Nozzle 
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distance from the nozzle exit plane. Figure 2.7 presents the IDA and ODA results compared 
with the measured data. 

The ODA methodology tends to underpredict the measurements in the near field of the 
plume. The flowfield that was used for these calculations started at the exit plane. Subsequent 
calculations that included the nozzle flowfield, as well as the plume, better predict the observed 
trends in the measurements. These calculations reproduce the test data beyond 2 feet and are 
20% low at 1 foot. The IDA results generally overpredict the data (up to 25-30%). It is 
possible that the boundary conditions that are used at the plume boundary are not appropriate 
for this application (although they are very good for radiation to the internal portions of nozzles). 
Further, research is required to investigate the potential effects on the IDA results due to 
boundary condition treatment. In view of the limited amount of validation that was performed 
with the ID ARAD code as compared to the SIRRM and REMCAR codes, the results are 
encouraging. 

2.3 Other Solution Methods for the RTE, Including Two-Flux Models 

Although the method of spherical harmonics discussed in previous section appeared to 
be an attractive approach to predicting radiation heating from SRM plumes, it had not been 
applied to this problem prior to this study. Therefore, several radiation analysis codes from the 
literature were considered for use in this study. The SIRRM-II code (Ref. 2.6) contains an 
extensive data base for gaseous narrow band models and particle radiation, and two-flux and six- 
flux radiation models. However, the SIRRM flux models are written for fore, aft, and side-on 
radiation analysis only, so they are of little direct value for base heating analysis. The 
REMCAR code (Ref. 2.15) is the reverse Monte Carlo code written by REMTECH; it is very 
general and useful, but it is also slow because of the extensive calculations required. The 
GASRAD code (Ref. 2.16) describes gaseous emission and absorption from H 2 0, CO z , CO, 
and soot for axisymmetric or three-dimensional flowfield input by integrating along multiple 
lines of sight. Among the other solution methods for the RTE considered was the method of 
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discrete coordinates (Ref. 2.19). This method has been used extensively and successfully for 
describing furnace operation. Furthermore, the first-order discrete coordinate solution is the 
two-flux model which forms the basis of the JANNAF SIRRM code (Ref. 2.6) and has been 
used in early work for SRM plume heating analysis. Therefore, a parallel study was made to 
determine the utility of the two-flux as an alternative analysis for SRM plume radiation 
evaluation. This two-flux study also provided a convenient tool for utilizing the extensive 
radiation property data base which already exists in the SIRRM code. The NOZZRAD code 
(Ref. 2.27) utilizes a two-flux model to describe emitting/absorbing/scattering media for an 
axisymmetric flowfield input. Gas and particle radiation are treated independently, not 
simultaneously in the NOZZRAD code. Both the NOZZRAD and GASRAD codes were used 
for sooty plumes. After the GASRAD or NOZZRAD code is used to establish the directional 
emissivities at the plume boundaries, the RAVFAC code (Ref. 2.23) is used to calculate 
radiation to points outside the plume. When RAVFAC is used to determine base heating, the 
view factors predicted with this code account for shading of vehicle components along the 
various lines-of-sight. Detailed descriptions of all of the codes mentioned in this paragraph are 
described in the cited references, except for the NOZZRAD code which is described herein. 

For emitting/absorbing media, integrations along lines-of-sight can be performed to 
predict radiation, as is done in the GASRAD code. If the media also scatters the radiation, the 
entire radiating volume must be considered at one time. If the volume consists of plane layers, 
each of which have constant properties, the radiative transport becomes essentially one- 
dimensional and the two-flux radiation analysis applies. Since the two-flux model resembles the 
gas only analysis, the same type of one-dimensional beam analysis can be applied to 
emitting/absorbing/scattering media if the following assumption is made. If the radiation field 
is assumed to be represented by a series of plane uniform layers which overlap and vary along 
each line-of-sight, multiple two-flux analyses can be performed to evaluate the local directional 
emission from the radiating volume. This procedure was used in the Aeronutronic work (Ref. 
2.28) and in the SIRRM code for field-of-view calculations. In SECA’s launch stand design 
studies (Ref. 2.27), the two-flux model analysis for slabs of varying temperature and particle 
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properties was updated and issued as the NOZZRAD code. In this study, the NOZZRAD code 
was extended to use an axisymmetric plume or motor analysis as input and perform one- 
dimensional beam analyses along selected lines-of-sight to provide directional emissivities from 
the plume. The resulting NOZZRAD predictions could then be used with RAVFAC to provide 
plume heating analyses. Thus, an analysis analogous to the GASRAD/RAVFAC predictions for 
gas plumes can now be performed for SRM plumes. The NOZZRAD analysis is developed as 
follows. 

2.3.1 The Two-Flux Model for Particle Flows 

The equation of radiative transfer along a line of sight: 

dl x {s,n,4>}/ds = -(a,+<T t ) I x {s,n,<f>} + a. U + (<x./4t) J 0 2t J o’ 

I x{s, /*.<£} p sin ©' dO'd <f>' (23) 

A beam of light which traverses an inhomogeneous medium is attenuated, a process 
called extinction, both by scattering of the light into other directions and by absorption. At a 
distance R from the scattering particle the scattered light has the character of a spherical wave. 
The direction of the scattered light is characterized by the angle 6 with the direction of the 
incident beam, and by the azimuthal angle <f>. The scattered intensity may be written as: 

I = I o P{0,0}/k 2 R 2 (24) 

where k = 2x/X is the wave number, F is the scattering function. If the total energy scattered 
is equated to the energy incident on an effective area <r„ it follows that: 

<r, = (1/k 2 ) J P{0,<£}du (25) 
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where do> = sin0d0d<£ and the integration is performed over all directions. The energy absorbed 
by the particle is set equal to the energy incident on the area <r„ and the total energy is set equal 
to the energy removed by the area a c . Therefore: 

<r e = a, + o; (26) 

At a given wavelength, the scattering and absorption cross sections of a spherical, 
homogeneous particle depend on only two parameters: the ratio of particle circumference to 
wavelength X = 2irr p /X, and the complex index of refraction m = n, - in 2 . For spherical 
particles of arbitrary size, all three of these cross sections can be determined by Mie theory. 
Since the scattering function is also determined by the Mie theory, the fraction of light scattered 
in a backward direction b is also determined. Tabulated values of n, and n 2 as a function of 
particle temperature T p , r p , and X are provided to a Mie code to yield <j„ and Eq. (29a). In 
fact, fractions of scattered radiation in any angle can be determined; for a six-flux radiation 
calculation fractions in the backward, forward, and sideways directions are so determined. 

Average values of <x, and a, over a particle size distribution are used, where o- = E N p 
<r i p /N, and the summation is on the particle size classes p. N t is the total number density of 
particles. Let ds = dz/fi, where n = cos 9. 

H dl x /dz = -N, (<r a +<rj I x + N, a, 1^, 

+ (N, cJAt) f 0 2r I x P d/x d 4 (27) 

Eliminate the phase function, P, using the "one-dimensional beam" approximation. 

dl x + /dz = -N t (a.+baJI x + + + bN t a.I x (28) 

-dlx/dz = -N t (<r,+bffJI x + N.a.I^ + bN t a,I x + (29) 
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where b = E [N p <r ap b p ]/(o, N,) and (29a) 

U = E [N p <r i p Ixb{T p }]/(<r a N.) 

Introduce the Optical Depth, r: 

dIx + /N t (<r a +cr,)dz = -I x + + [(l-b)a,/((T a +crJ]I x + (30) 

+ [a a /(a a +aJ]I Xb + [bay(a a +aj]l x 
b = fraction back scattered radiation = 1-f 
f = fraction forward scattered radiation 
dr = N t (a a +<r s )dz, the differential optical depth 
where a x = [a/fo+aj], the albedo 
1- a x = K/K+aJ] 

Wavelength dependence of f & b will not be indicated. 

+dl x + /dr = -I x + +(l-b)a x l x + +(l-a x )l xb +ba x l x (31) 

+dl x + /dr = -I x + +a x (fl x + +bl x )+(l-a x )l xb (32) 

-dI x -/N t (<T a +<r s )dz = -I x +(l-b)a s I x /( ( j a +a s ) (33) 

+bcr s I x + /(<r a +a 1 ) + <r a I Xb /(a a +cO 
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-dlx/dr = -I x +a x (fI x -+bI x + ) + (l-a x )I Xb (34) 

Assuming all intensities and optical properties are monochromatic X will be omitted from 
here on. 

+dl + /dr = -I + + a(fl + +bl) + (l-a)I b (35) 

+dI7dr = +r - a(fl'+bl + ) - (l-a)I b (36) 

Let m = cr,/(a,+a,)=(l-a) or a = 1-m 

k = [m(2b(l-m)+m)]° 5 (37) 

k = m[(2b(l-m)/m) + l] 05 
(k/m) I 2 = [(2b/m)(l-m) + l] 

dl + /dr = -0.5[(k 2 /m)+m]I + 0.5[(k 2 /m)-m]L + ml b (38) 

dl'/dr = +0.5[(k 2 /m)+m]J - 0.5[(k 2 /m)-m]I + - ml„ (39) 

The solution of the equation of transfer for a single isothermal slab is: 

I + = 0.5A[(k/m) + l]e kT + 0.5C[(k/m)-l]e kT + I b (40) 

T = 0.5A[(k/m)-l]e kr + 0.5C[(k/m)+l]e^ + I b (41) 

This solution for several isothermal slabs is obtained by imposing the interface boundary 
conditions: 

I + m-i {r,} = I + m( T i} & r*,{r,} = T m {rd (42) 
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If the interface position, i, and nearest slab, m-1, are numbered the same, say p, and if the next 
slab out is numbered q, these conditions become: 

I\ M = & r p {r p } = r,{r„} (43) 

Surface boundary conditions: 

i+ o{t 0 } = e 0 I b {r 0 } + r 0 I'JrJ (44) 

= e a I b {r n } + r n I + n {r n } (45) 

where the emissivity, e, and reflectivity, r, are for the environments of the slabs. 

With these boundary and interface conditions, the solution for several isothermal slabs is given 
by: 

I + 0 = O.SA.Kk/mHl], + 0.5C,[(k/m)-l], + I bl (46) 

I + „{r p } = 0.5A p [(k/m) + l]pe kp Tp + 0.5Cp[(k/m)-l]pe** + I* (47) 

I + q {r p } = 0.5A q [(k/m) + l] q e kqTp + 0.5C q [(k/m)-l]^^ + I,, (48) 

Ip{r p } = 0.5 A p [(k/ m)- l] p e' kp ^ + 0.5C p [(k/m) + l] p e kpTp + I* (49) 

r q {r p } = 0.5A q [(k/m)-l] q e' kqrp + 0.5C q [(k/m) + l] q e k< ^ + I,, (50) 

I n = 0.5A n [(k/m)-l] n e kn ” > + 0.5Q.[(k/m)+l] 1 jB>— + I to (51) 


2-78 



SECA-FR-94-18 


If the I + ’s and I’s at the interfaces are eliminated, there are sufficient equations to evaluate the 
A’s and C’s. The solution of these equations is provided by the NOZZRAD code. 

2.3.2 The NOZZRAD Code 

NOZZRAD was developed by SECA to predict radiation heat transfer to a point on a 
nozzle wall or plume boundary. NOZZRAD has also been set up to be used in conjunction with 
the RAVFAC (Ref. 2.23) code for prediction of radiation heat rates to surfaces outside the 
flowfield boundaries. 

Flowfield data for NOZZRAD is supplied in an input file of the same format used in 
SIRRM (Ref. 2. 16) and must be axisymmetric. The nozzle wall or plume boundary is assumed 
to be the outer boundary of this flowfield map. Lines-of-sight (LOS), from a specified point on 
the nozzle wall or plume boundary, are drawn across the flowfield at evenly spaced angular 
intervals. Flowfield properties are obtained at specified distances along each specific LOS 
creating one-dimensional slabs from which specific intensity is calculated from the two-flux 
method described in Section 2.31. The specific intensities for each LOS are appropriately 
integrated to calculate the total and average radiation intensity to the specific point. 

NOZZRAD has the capability to calculate radiation from either A1 2 0 3 or carbon/soot 
particles and the gaseous species of H 2 0 and C0 2 . Particle and gas calculations are treated 
separately. Two options for the gaseous radiation calculations are included. The first option 
treats the gas as one isothermal, homogeneous slab by averaging the points along the LOS. The 
second option treats the composition and temperature across the LOS as a summation of 
isothermal slabs. 

A1 2 0 3 particle optical properties are read from the files SIGSCL01.DAT, 
SIGACL01.DAT, and BETACL01.DAT. Carbon/soot particle optical properties are read from 
the files SIGSCL02.DAT, SIGACL02.DAT, and BETACL02.DAT. These files contain 
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scattering and absorption cross-sections and back scatter efficiencies assuming a spherical shape 
as a function of 10 particle sizes, 5 particle temperatures, and 37 wavelengths. More 
information on these files can be found in Section 2.1. 

Flowfield properties are read into NOZZRAD via a file named SIRRM.DAT. This file 
is a standard SIRRM flowfield map which is generated by the Standard Plume Flowfield Code 
(SPF/2) (Ref. 2.29) with five gaseous species in the following order: H 2 0, C0 2 , HC1, CO, OH. 
More information on this flowfield data file can be found in Section 2.1. 

User inputs are read into NOZZRAD via a file named NOZZRAD.INP. This file 
controls all of the user options available for NOZZRAD. These options include: 

1) Case or NOZZRAD Run Description 

2) Type of Particle or Gas Radiation Calculation 

3) Number of Angular Intervals 

4) Thickness of Slabs along a LOS 

5) Fields of View 

6) Angles of Orientation . 

An example NOZZRAD.INP is given in Table 2. 16. The format of NOZZRAD.INP 
along with the descriptions of each input variable is given in Table 2. 17. The source code for 
NOZZRAD, along with the required data files and sample input files, are contained on the MS- 
DOS disk, RAD7. 
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Table 2.16 Example of Input File - NOZZRAD.INP 


EXAMPLE OF NOZZRAD USER INPUT 

lOPT-flow typo option, 1 » AL203 particle,2 ~ Soot particle, 3 =gas 4«gas avg 

1 '-i ; '• 

M -number of LOS intervals 
5 

ODD -LOS INTERVAL {ft) 

.0500 :: ; 

K-wall boundary: 1 -wall perpendicular to axis, 2 -radial wall 

2 ■ : : 

Number of J-POINTs to evaluate and number of the J-POINT 
9 

4,7, 1 2, 1 5,4,5,6,14,9 
FOV -Field of View (deg.) 

0, l 0.,0. / 0.,0. ( 0.,0.,0.,0.,0.,0. r 0.,0.,0.,0.,0.,0..0.,0. 

Orientation option: 1 = user defined, 0 * perpendicular to wall or plume boundary 

■ ; : i h - v: : ■ ; : ■ : : : i : ^ ; : : : : -; : v^ 

SLP -User defined orientation ( + deg) 

0 M 0.>0 ,0 , 240 r ,30 ,, 60 ,30 , 50., 




SECA-FR-94-18 



Description 




Case or Run Description 


Dummy variable to describe next input 


Puticjc/Gu Radiation Calculation Option 


0 = Particle Radiation Calculations for ALjOj 

1 = Particle Radiation Calculations for Carbon/Soot 

3 - H a 0 &. CO, Gaseous Radiation Calculations -One Slab, Average Properties 

4 = HjO & COj Gaseous Radiation Calculations -Multiple Slab 





Dummy variable to describe next input 


1/2 the number of angular subintervals to be used in the composite numerical integration scheme 
(Simpson’s Rule). Radiation intensities for (2M + 1) 2 lines of sight will be calculated. (MAX = 100 ) 


Dummy variable to describe next input. 


Distance between points along the line of sight where flowfield properties are "looked up". For particle 
radiation, this distance can also be thought of as the thickness for the slabs in the radiation calculations. 


Dummy variable to describe next input. 


1 = Inlet Wall (Nozzles) or Exit Plane (Plumes) 

2 - Nozzle Radial Wall or Plume Radial Boundary 

(See Figure 1.) 




Dummy variable to describe next input. 


Number of wall points for which radiation calculations will be made. 


Identity of each wall point. These are integer values corresponding to the boundary points in the flow fie Id 

data file, SIRRM.DAT 




Dummy variable to describe next input 


Field of View 

Full field of view seen by each wall point Corresponds to the same order of the J-POINT wall 

identification. 

(See Figure 1) 






Dummy variable to describe next input 


User Defined Orientation Option 


0 = Field of View is oriented normal to the specified wall. 

1 = Field of View is oriented as defined by user. 


Dummy variable to describe next input 


User Defined Orientation Angles 

Orientation angles for each specified wall point if ISLP= 1. Corresponds to the same order of the J- 
POINT wall identification. 

(See Figure 1) 
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The RAVFAC (Ref. 2.23) code can be used in conjunction with NOZZRAD for the 
calculation of radiative heat flux to surfaces outside a rocket plume boundary. RAVFAC is 
essentially a view factor code which represents a rocket plume as a geometrical surface with 
specific emissive intensities. With the proper input, NOZZRAD will generate all RAVFAC input 
files associated with the plume surface. These files can be run with the original RAVFAC code 
which assumes diffuse surface radiation or with a modified version of the RAVFAC code which 
accounts for variations in directional intensities. The modified RAVFAC code reads and 
properly evaluates directional intensity data stored on an additional file called SPECINT.DAT 
which is generated by NOZZRAD. The diffuse surface radiation assumption should predict 
conservative answers in relatively short NOZZRAD run times. Directional considerations should 
be more accurate; but, they require more time for NOZZRAD to calculate. Table 2. 18 provides 
the user with a guide for setting up the appropriate geometrical inputs to NOZZRAD. INP to 
generate inputs for the original RAVFAC code (diffuse emission). Table 2. 19 provides the user 
with a guide for setting up NOZZRAD. INP to generate inputs for the modified RAVFAC code 
(directional emission). 

The output file named RWALL.DAT provides the user with the total radiation heat flux 
and average radiation intensities as calculated for each user specified boundary point. The total 
radiation heat flux prediction accounts for the normal surface to LOS angle difference to provide 
predictions which are comparable to radiometer values. Intensities are averaged over the user 
specified field of view. An example of RWALL.DAT is given in Table 2.20. 
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Table 2.18 NOZZRAD Input for use with RAVFAC (diffuse emission) 


Variable 

Description 

K 

USE 2 

NL 

The number of boundary points (J -Points) in the axial direction needed for the accurate 
description of the plume surface and emissivity. Generally the user should use the total 
number of axial points described in SIRRM.DAT. 

J 

Identity of axial boundary point. For the proper generation of the plume surface description 

these must be in consecutive order. 

FOV 

USE 0 for each J 

ISLP 

USE 0 

(LOS normal to plume surface) 


Table 2.19 NOZZRAD Input for use with modified RAVFAC (directional emission) 


Variable 

Description 

M 

User should chose this variable based on the angular increment which is desired between 
directional intensities. The suggested value of 1 1 will provide a maximum angular increment 

of 18°. 

K 

USE 2 

NL 

Number of boundary points (J-Points) in the axial direction needed for the accurate 
description of the plume surface and emissivity. Generally the user should use the total 
number of axial points described in SIRRM.DAT. 

J 

Identity of axial boundary point. For the proper generation of the plume surface file, these 

must be in consecutive order. 

FOV 

USE 180 for each J 

ISLP 

USE 0 

(LOS normal to plume surface) 
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Table 2.20 Example of the Output File RWALL.DAT 


♦A****************************************************************************** 

SECA"S NOZZRAD 


AL203 PARTICLE RADIATION CALCULATIONS 
H- INTERVAL « 5 

■ K-WALL * 2 : ' 

J -POINTS = 4, 7, 12, 15, 4 

5, 6, 14, 9 , 

LOS INCREMENT * 0.5000E-01 (ft) 
******************************************** 

EXAMPLE OF NOZZRAD OUTPUT 


RADIATION HEAT RATE TO WALL or PLUME BOUNDARY 


(BTU/f t2/S) 

X (ft) R (ft) RADIATION HEAT RATE AVG. INTENSITY FOV (deg) S (deg) 


0.9537E+00 0.1353E+01 0.70772E+02 0.22893E+02 .00 .00 


The output files PLUME.DAT, ICS.DAT, and RAD.DAT are standard input for the 
radiation view factor code RAVFAC (Ref. 2.23). The output file named SPECINT.DAT is a 
data file with the specific intensities of every calculated line of sight. This file provides the 
directional intensities for the modified version of RAVFAC described previously. 

2.3.3 Results of the NOZZRAD Analysis 

Results of the MNASA Test Measurements along with comparisons to SIRRM-II and 
Monte-Carlo based predictions were used to validate the NOZZRAD methodology. The plume 
flowfield which was used for the SIRRM-II, Monte-Carlo, and NOZZRAD predictions was 
generated with SPF/2 and is discussed in detail in Ref. 2.29. The MNASA test setup is 
illustrated in Fig. 2.8. Radiometer instrument numbers 1-9 are small field of view instruments. 
Radiometers 17-21 are full field of view instruments. Table 2.21 shows the experimental results 
for radiometers 1-9 along with the predictions made from SIRRM-II, Monte-Carlo and 
NOZZRAD. NOZZRAD was run in these cases for A1 2 0 3 particle radiation with one line- of- 


2-85 





SEC A-FR-94- 1 8 



SYMBOLS: 

Narrow-view radiometer 
• Wide-Angle radiometer 


Fig. 2.8 Radiometer Orientations for the MNASA Tests 
(Ref. REMTECH RTN 213-18) 
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sight. Both SIRRM-II and NOZZRAD intensity values were multiplied by the cosine of the 
angle between the surface of the plume and the radiometer to provide the correct view factor. 
The MONTE-CARLO predictions were taken from (Ref. 2.30). With the exception of 
Radiometer #5 the results obtained from NOZZRAD were reasonable. 

Table 2.21 Comparison of MNASA Test Measurements for Small Field of View Radiometers 
with MONTE-CARLO, SIRRM-II, and NOZZRAD Predictions 


INSTRUMENT 

NUMBER 

TEST 

(Btu/ft J /s) 

MONTE-CARLO 

(Btu/ftVs) 

SIRRM-II 

(Btu/ftVs) 

NOZZRAD 
AL^Oj Particle 
(Btu/ftVs) 

3 

6 

8 

11 

1 

55-56 

55-57 

58-64 

- 

67.1 

64.8 

66.7 

2 

44 

46 

50 

43-48 

52.8 

40.1 

45.5 

3 

38-43 

46-48 

50-54 

- 

57.4 

47.2 

56.9 

4 

10-29 

36-46 

38-47 

29-30 

58.1 

45.7 

55.2 

5 

66-67 

51-56 

59-60 

- 

99.5 

60.8 

87.4 

6 

50-53 

47-49 

56-57 

- 

54.2 

53.7 

50.0 

7 

43 

38 

38 

- 

48.6 

58.2 

54.2 

8 

35 

41 

43 

- 

60.0 

46.6 

57.8 

9 

38-43 

34-39 

37-42 

- 

54.5 

35.3 

49.2 


Since gaseous and particle radiation is treated separately in NOZZRAD, similar 
predictions were made with SIRRM-II where gaseous and particle radiation was separated by 
altering the flowfield input. Gas partial pressures were set to zero for A1 2 0 3 particle radiation 
predictions and particle number densities were set to zero for gaseous radiation calculations. 
These calculations were made for a single line-of-sight for each of the Radiometers numbers 1-4 
of the MNASA tests. Results for the gaseous radiation predictions are shown in Table 2.22. 
Results for the particle radiation calculations are shown in Table 2.23. Good agreement is 
shown for the gaseous results but the particle radiation results show relatively large differences 
even after the regular NOZZRAD properties were replaced by the same properties used in 
SIRRM-II. Further investigation revealed that these differences are due to the differences in the 
optical property interpolation schemes used by the two codes. Figure 2.5 showed the spectral 
emissive power of a 4.5 ft thick homogeneous slab of A1 2 0 3 particles with radius of 3 microns 
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and 100,000 cm' 3 in number density. Particle temperatures were varied between the tabulated 
values of 3000°K and 2320°K. At the tabulated values there is excellent agreement; but, 
NOZZRAD predicts higher values at temperatures in between the tabulated temperatures. It is 
believed that similar differences will be present as a function of particle size but specific 
calculations in this regard have not been made. No inferences should be made from Fig. 2.5 
concerning the conservative predictions made by NOZZRAD. Non-homogeneous slab 
configurations can be combined to produce NOZZRAD predictions which are lower than 
SIRRM-II predictions as evident from Table 2.23. The reason that the particles only SIRRM 
solution is higher than the particles, plus gas, is shown in Fig. 2.6 and is due to the assumed 
gas particle radiation interaction coded in SIRRM. This effect is very large for this case. More 
consideration of the gas/particle interaction and of the interpolation method used to obtain optical 
properties of the particles is evidently needed to reconcile these large differences in these 
prediction methods. 


Table 2.22 Comparison of SIRRM-II and NOZZRAD for Gaseous Radiation Predictions 


INSTRUMENT 

ft 

SIRRM-II 

(Btu/fF/s) 

NOZZRAD 
Multiple Slab 
(Btu/ft 2 /s) 

NOZZRAD 
One Averaged Slab 
(Btu/ff/s) 

1 

10.1 

10.3 

8.6 

2 

12.98 

13.0 

10.6 

3 

19.7 

17.9 

13.2 

4 

16.7 

16.9 

11.8 


Table 2.23 Comparison of SIRRM-II and NOZZRAD for ALj 0 3 Particle Radiation Predictions 


INSTRUMENT 

SIRRM-II 

NOZZRAD 

NOZZRAD 

ft 

(Btu/ft 2 /s) 

(Btu/fE/s 

SIRRM Properties 

1 

69.2 

66.7 

70.8 

2 

64.2 

45.5 

48.2 

3 

73.8 

56.9 

60.3 

4 

71.8 

55.2 

58.5 
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In order to predict radiation heat flux to surfaces outside the plume boundary with large 
fields of view, NOZZRAD should be used in conjunction with the RAVFAC code. The results 
of this methodology as applied to the MNASA test data are shown in Table 2.24. NOZZRAD 
was run according to the guides provided in Tables 2.18 and 2.19 for A1 2 0 3 particle radiation. 
As expected the diffuse assumption provided conservative answers. Radiometers 17 and 18 were 
much more sensitive to directional considerations due to their orientation with respect to the 
plume surface (see Fig. 2.8). The fields of view for these radiometers see intensities with 
angular directions far from the normal of the plume surface. 

Table 2.24 Comparison of MNASA Test Measurements for Full Field of View Radiometers 


with RAVFAC/NOZZRAD methodology 


INSTRUMENT 

NUMBER 

TEST 

(Btu/fF/s) 

RAVFAC 

DIFFUSE 

(Btu/ft 2 /s) 

RAVFAC 

DIRECTIONAL 

(Btu/fP/s) 

3 

6 

8 

11 

17 

2.80 

2.66 

2.75 

2.83 

4.63 

2.77 

18 

3.11 

3.00 

3.07 

- 

4.25 

2.99 

19 

5.47 

5.27 

5.39 

- 

6.22 

5.21 

20 

7.12 

6.99 

7.21 

- 

8.46 

7.30 

21 

9.73 

9.30 

9.54 

- 

11.15 

9.99 


Other general comments on the use of NOZZRAD: 

1) The number of slabs for a particular LOS is limited to 200. 

2) Numerical integration errors for the total radiation heat 

flux calculations may arise if the number of angular increments 
which is chosen is too low. This problem occurs to a greater 
extent as the field of view angle becomes smaller. For field of view 
angles < 30° it is suggested that the user assume diffuse radiation 
and chose an FOV of 0°. 


NOZZRAD has been validated with some comparisons to SIRRM (Ref. 2.6) and general 
isothermal slab solutions as well as some experimental data, however; until NOZZRAD has been 
used more extensively caution should be exercised with the use of the results. 
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2.3.3 Radiation from Sooty Plumes 

Current NASA interest in hybrid motors and RP-1 fuel has prompted including soot 
radiation in this study. Soot radiation can be predicted with either NOZZRAD or the GASRAD 
code (Ref. 2.16). Soot particles in plumes are believed to be small enough that they do not 
scatter radiation. This assumption is made in the GASRAD program. However, the required 
mole fraction of soot is a difficult variable to evaluate. Not only must the chemistry of the 
sooting combustion be described, but the molecular weight of the soot must be specified. Such 
predictions are not within the scope of the current investigation; therefore, soot was 
approximated as a specified fraction of the carbon in the fuel with the thermodynamic properties 
of graphite. 

Preliminary analysis of radiation heating rates has been performed on a hybrid motor 
using SPF/2 predicted plume with 2% carbon. Radiation heating rates were calculated using 
NOZZRAD and GASRAD for a full field of view at various points downstream along the plume 
boundary. Figure 2.9 illustrates the positions and calculated heat rates from NOZZRAD using 
a particle size of 0. 1 micrometer. The radiation heating rates calculated by GASRAD for the 
same points were negligible. When GASRAD was run without the cool outer layers of the 
plume, the radiation heat rates which were calculated were comparable to NOZZRAD. 
Apparently, the assumed carbon content in the low temperature shear layer absorbed most of the 
radiation from the high temperature portions of the plume. The discrepancies in the answers 
from NOZZRAD and GASRAD can be attributed to particle size. GASRAD assumes the carbon 
particles are so small that no scattering effects are present. Since radiation heating rates which 
are shown in Fig. 2.9 are reasonable in comparison with data from similar engines, this initial 
investigation indicates that scattering for the soot particles should be considered for the assumed 
carbon distribution used in the SPF/2 plume flowfield prediction. However, the real problem 
is to accurately predict the soot concentrations and particle size distributions. 
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Fig. 2.9 Plume Radiation from Soot 
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3.0 TWO-PHASE FLOW MODELING FOR 
SOLID ROCKET MOTOR RADIATION PREDICTIONS 

The spatial characterization of gas and particulate properties of solid rocket motor nozzle 
and exhaust plume properties is more important than the radiation models which are used to 
determine the radiation fluxes that are emitted from the flowfields. Even if the radiative model 
exactly models all the radiation processes of the gas and particulates, it is impossible to perform 
a radiation prediction if the flowfield is improperly characterized. Thus an important part of the 
investigation of new techniques for solid rocket motor radiation predictions was the investigation 
of the adequacy and accuracy of the existing models which are available to predict solid rocket 
motor flowfields. This section of the report describes the results of the evaluation of solid 
rocket motor flowfield models. In addition to the actual flowfields models, submodels such as 
particle size models (Section 3.3), and soot (Section 2.3.3) were investigated relative to the 
importance of the submodels used by the flowfield codes in predicting radiation. 

3.1 Conventional Solid Motor Flowfield Prediction Methodology 

The most commonly utilized model for calculating solid rocket motor flowfields for low 
altitude solid rocket motors is the JANNAF sponsored Standard Plume Flowfield Model. The 
older versions of the SPF code (Ref. 3.1) (SPF1 and SPF2) required that the combustion 
chamber-nozzle flowfield be calculated with another code and passed to the SPF code in the 
form of exit plane distributions of gas and particle flow properties. A typical code used to 
supply exit plane properties to SPF is the RAMP2 code (Ref. 3.2). 

The RAMP2 code has been continuously improved under NASA funding since the mid 
70’s. This code was originally developed to support the Space Shuttle design studies for the low 
to mid altitude flight regions. In the early 80’ s the capabilities of the code were extended so that 
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vacuum plumes could be treated for in orbit spacecraft design applications. The code calculates 
an inviscid flowfield and is limited to flows which contain no imbedded subsonic regions. 

Rocket exhaust flowfields are very complicated and are governed by many phenomena. 
Many simplifying assumptions are made to enable one to compute exhaust flows. However, 
many of these simplifying assumptions can compromise and invalidate the results, depending on 
the application for which the flowfield is intended. Numerous inviscid codes are available that 
treat many of the governing phenomena, but no single code is available that treats reacting 
single- and multi-phase flows including boundary-layer effects as an integral part of the solution. 
Thus, previously it was necessary to use a multitude of codes to treat inviscid nozzle/plume flow 
in detail. It is therefore desirable from both computational and economic standpoints to have 
a single code that can treat all the dominant phenomena in a rocket nozzle/plume flowfield. 
Additionally, it is possible to perform calculations which may range from the most simple (as 
for preliminary design studies) to the most complex as required for final design. 

The basic RAMP2 code employs modular construction and has the following capabilities: 
(1) Two-phase with a two-phase transonic solution, (2) Two-phase, reacting gas (chemical 
equilibrium, reaction kinetics), supersonic inviscid nozzle/plume solution, and is (3) Operational 
for inviscid solutions at both high and low altitudes, (4) Direct interface with the JANNAF SPF 
code, (5) Shock capturing finite difference numerical operator, (6) Two-phase, 
equilibrium/frozen, boundary-layer analysis, (7) Variable oxidizer-to-fuel ratio transonic 
solution, (8) Improved two-phase transonic solution, (9) Two-phase real gas semi-empirical 
nozzle boundary layer expansion, (10) Continuum limit criteria, and (11) Sudden freeze free 
molecular calculation beyond the continuum limit. 

Most of the above capabilities already exist in other computer codes. These codes were 
incorporated into the RAMP code to enhance its usefulness. 

The three programs which make up the RAMP2 code (TRAN72-Ref. 3.3), RAMP2F, 
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and (BLIMP J-Ref. 3.4) have been modified so as to interact as if they were a single code even 
though they are executed separately due to computer storage restrictions. 

In general, in order to solve a high altitude plume the following steps are required. First, 
the TRAN72 program input data is prepared and executed to generate a data file describing the 
thermodynamic characteristics of the post-combustion gases. Next, the RAMP2F flowfield data 
are prepared and the nozzle flowfield is solved using the TRAN72 program data file as input. 
Then, in order to adequately describe the nozzle boundary layer, the BLIMPJ code is executed 
using an input data file and the flowfield file generated by the RAMP2F nozzle solution. Finally 
the exhaust plume is calculated by using the nozzle solution and boundary layer solution to 
generate an exit plane start line that is used to initiate the plume solution. Thus, the generation 
of a high altitude plume can require up to four different executions of programs (TRAN72, 
RAMP2F, BLIMPJ, and RAMP2F) for the specification of the most detailed and accurate 
results. Physical input data are required only for the TRAN72 and first RAMP2F execution. 
All data required for the BLIMPJ code and second RAMP2F execution are generated internal 
to the program and/or communicated via data tapes or temporary files. Depending on the 
application, the problem, or the level of sophistication required in the plume results, it may not 
be necessary to run the TRAN72 or BLIMPJ codes. It is possible that a single RAMP2F 
calculation may be adequate, such as in the case for a low altitude plume, which is what was 
done in this study to support the ASRM flowfield modeling. For low altitude cases the RAMP2 
code was used to generate the exit plane start line which is used by SPF/2 to initiate the plume 
solution. 

The Joint-Army-Navy-NASA-Air Force (JANNAF) Standard Plume Flowfield (SPF) 
Model is a modular computer program which has been under development several years by 
Science Applications International Corporation (SAIC) of Wayne, PA. The development of this 
program has been sponsored by the U. S. Army Missile Command, (AMICOM) at Huntsville, 
AL, NASA at Langley Research Center, VA, and Arnold Engineering Development Center 
(AEDC) at Tullahoma, Tennessee. The program has undergone three stages (SPF/1 , SPF/2, and 
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SPF/3) in its development. 

The version of the SPF Program that was used in this study to investigate all flowfields 
is the SPF/2 version. In SPF/2 the input is supplied both from the user and from the data bank 
component. The data bank contains JANNAF thermodynamic data from about 95 chemical 
species and chemical reaction data from about 107 chemical reactions. This input data then goes 
to the processor component, PRC2 of the program. The output from PRC2 serves as input to 
the inviscid component, SCP2 (provided a 2-D startline is desired). The output from SCP2 is 
then used as input to the plume mixing layer component, BOT2. If a 1-D startline is desired, 
no inviscid calculation will be obtained, and the output from PRC2 will go directly to BOT2. 

The SPF/2 Program has the capability of treating the following six chemical systems: 
1) H/O, 2) C/H/O, 3) C/H/O/Cl, 4) C/H/O/Cl/F, 5) H/O/B, and 6) H/O/B/Cl/F. In addition, 
another system may be used in which the user specifies the chemical species. 

For a 2-D input across the exit plane (or separation plane) the input was obtained from 
the output of the RAMP2 program with a distribution of the gas temperature, pressure, axial and 
radial velocity; particle density, velocity and temperature at each radial point on the startline. 
The chemical species are frozen across the exit plane and for the entire length of the inviscid 
plume. 


The plume flowfield generated by SPF/2 is calculated by SCP2 for the internal inviscid 
core (hyperbolic solution). The outer annulus or plume mixing layer (parabolic solution) is 
computed by BOT2. 

The SPF/2 program is used primarily at the lower altitudes where the Mach disc is an 
important contribution to the overall base radiation flux and where mixing and afterburning along 
the plume boundary play an important role in the base environment. 


3-4 



SECA-FR-94-18 


RAMP2 and SPF2 codes have been extensively used to perform solid rocket motor 
flowfield predictions that were subsequently used by radiation codes to predict radiation fluxes 
to vehicle structure. The RAMP2 and SPF2 codes have been improved (Ref. 3.5) to the point 
that along with the NASA funded REMCAR (Ref. 3.6) radiation code accurate prediction of 
radiation loads to launch vehicles and missiles are now possible. These improvements in 
RAMP2 and SPF2 are referred to as the Cycle 2.0 methodology. 

3.2 Two-Phase Navier Stokes Flowfield Modeling 

Navier Stokes flow solvers have reached a level of maturity that potentially could result 
in two-phase flowfields which could be utilized to perform radiation predictions of launch vehicle 
plume induced radiation heating. Under a previous NASA funded study (Ref. 3.7), a particulate 
two-phase model was incorporated into an existing, gas only, Navier-Stokes Computation Fluid 
Dynamics code (CFD). The code which was used as the basis of the new code was the FDNS 
code (Refs. 3.8, 3.9 and 3.10). 

The FDNS code solves a set of nonlinear and coupled transport equations, the Navier- 
Stokes equations, energy equation, two-equation turbulence models and chemical species 
continuity equations in non-dimensional form. Finite difference approximations are employed 
to discretize the transport equations on non-staggered grid mesh systems. High-order (second- 
or third-order) upwind or central differencing schemes plus adaptive second-order and fourth- 
order dissipation terms are used to approximate the convective terms of the transport equations. 
Second-order central differencing schemes are used for the viscous and source terms of the 
governing equations. To insure positive numbers for some scaler quantities such as turbulence 
kinetic energy and species mass fractions, a first-order upwind scheme is employed for the 
convection process. A pressure based predictor/multi-corrector solution procedure is employed 
in the FDNS code to enhance velocity-pressure coupling and mass-conserved flowfield solutions 
at the end of each time step. This pressure based method is suitable for all speed flow 
computations. A time-centered Crank-Nicholson time-marching scheme is used for the temporal 
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discretization for time-accurate solutions. For steady-state flows, an implicit Euler time-marking 
scheme can be used for better computational efficiency. The selection of time marching scheme 
can be used for better computational efficiency. The selection of time- marching schemes is 
controlled through input data. 

In the current version of the FDNS code (Ref. 3.11), incompressible or compressible, 
standard or extended, k-e turbulence models with wall function or direct integration to the wall 
(low-Reynolds number turbulence model) options are included. Turbulence model options are 
selected through input data. Chemical kinetics and species thermodynamics data are required 
to be prepared in the input data file. 

For particulate two-phase flow simulations, a Lagrangian method using an implicit 
particle trajectory integration scheme is used. In the present version of FDNS, called FDNSEL, 
only steady state (not time-varying) solutions of two-phase flow is possible. This section of the 
report describes: theories that are incorporated in FDNS (3.2.1), the history and validation of 
FDNS two-phase flow version (3.2.2), FDNS input instructions and sample cases for two-phase 
nozzle analysis (3.2.3) and the possible influence of combustion chamber geometry on predicted 
radiation (3.2.4). 

3.2.1 FDNS Theories 

This subsection describes some of the basic theories that are incorporated into FDNS. 
More detailed descriptions of the theories and basic code description can be found in Ref. 3.11. 

Governing Equations . The gas-phase governing equations of the FDNS module are the 
Reynolds-averaged Navier-Stokes equations with the addition of particle drag forces and heat 
fluxes in the momentum equations and the energy equation, respectively. Due to the effect of 
large density differences between the particles and the surrounding gas, the drag force was 
considered to be the primary contribution to the inter-phase momentum exchange. The gas- 
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phase governing equations are written as: 


J-^dpq/at) = a[-pUq + n ef &s(d q/d^J/d*, + S q 


where q = 1, u, v, w, h, k, e and for the continuity, momentum, energy, turbulence model 
and chemical species transport equations respectively. And, the transformation parameters and 
effective viscosity, p eff , are given as: 


J = d(£,r?,f)/d(x, y,z) 
Uj = (Uj/JX^/dXj) 

Gy = mdxjidtydxjn 

Pcff = O + Pt)/tfq 


The source terms in the governing equations, S q , are given as: 


= /-* 


0 

~Px + v Kr(«,)J " (2/3) (n^vu) x + D 
“P, + v [p^(« / ) y ] - (2/3)(p^v«) y + D 

~Pz + v t/vW " ( 2/3 ) (/V™) z + D z 

DP/DT + h v + H p - u p D x - v p D y - w p D z 

p(P r ~ e) 

p (e!k)[( C, + C 3 P /e) P - C 2 e ) 


where Dx, Dy and Dz represent the drag forces and c takes on values between 1 and N (number 
of gas species). u p , v p and w p are the particle velocity components. H p is the rate of heat 
transfer per unit volume to the gas phase. hy stands for the viscous heat flux of the gas phase. 
P r stands for the turbulence kinetic energy production rate and is written as: 
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P r = (/i/p)[(3u i /3x i + du;/dXj) 2 /2 - 2(du k /dx k ) 2 /3] 


An equation of state, p = p/(RT/M w ), is used to close the above system of equations. Turbulent 
Schmidt and Prandtl numbers, a q , for the governing equations and other turbulence model 
constants given, are taken from Refs. 3.12, 3.13 and 3.14. 


Finite Rate Chemistry Model . For gas-phase chemical reaction modeling, a general system of 
chemical reactions is written in terms of the stoichiometric coefficients and v t] ') and the i-th 
chemical species name (MJ of the j-th reaction as 

E v 9 M, = E v { ; M/ 

i i 


The net rate of change in the molar concentration of species i due to reactions j, X^, is 
written as: 

Xij = Oy-^y) [K Q n(WM w ,r - JZfKfiaJ M Wi )' ij ] 

and the species production rate, coj , (in terms of mass fraction) is calculated by summing over 
all reactions. 

= Ki E ** 

j 


where 

M wi = molecular weight of species i 

ofj = mass fraction of species i 

p = fluid density 

Kg = forward rate of reaction j 

K bj = backward rate of reaction j = Kg/K^ 

Kej = equilibrium constant = (l/RT) E(,ij, '' ij) exp{E(f i V ij ' - f^jj)} 
f; = Gibbs free energy of species i 
K f = A T B exp{-E/RT} 

Finally, the species continuity equations are written as: 
p D t of ; - V[(p cfr /aJVaJ = Wj 

where o a (assumed to be 0.9) represents the Schmidt number for turbulent diffusion. Either a 
penalty function or an implicit integration is employed to ensure the basic element conservation 
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constraints at the end of every time marching step. This is a crucial requirement for the 
numerical stability and accuracy of a CFD combustion model. The penalty function calculation 
is accomplished by limiting the allowable changes in species concentrations, which are the 
solutions of the species continuity equations, for each time step such that the species mass 
fractions are well bounded within physical limits. The resulting limited changes are adjusted so 
that they are proportional to the species source terms. A similar chemistry approach and 
detailed turbulence submodels were reported previously (Ref. 3. 15). 

Particulate-Phase Equations . A Eulerian-Lagrangian particle tracking method is employed in 
FDNS to provide effects of momentum and energy exchanges between the gas phase and the 
particle phase. The particle trajectories are calculated using an efficient implicit time integration 
method for several groups of particle sizes by which the drag forces and heat fluxes are then 
coupled with the gas phase equations. The equations that constitute the particle trajectory and 
temperature history are written as: 


DV/Dt = (Uj - Vi)/t d 

Dhp/Dt = Cp, (T aw - T p )/t„ - 6 ad T p 4 /(p p d p ) 


where Uj = Gas Velocity 
Vi = Particle Velocity 

t d = Particle Dynamic Relaxation Time = 4 p p dp/(3 C d p c | U s - V s J ) 
h p = Particle Enthalpy 

= Particle Heat Capacity 
T p = Particle Temperature 
T, w = Gas Recovery Temperature 

t H = Particle Thermal-Equilibrium Time = (p p d p )/[12 Nu /*/(Pr d p )] 
a = Stefan-Boltzmann Constant = 4.76E-13 BTU/FT 2 -SR 
e = Particle Emissivity = 0.20 — 0.31 
f = Radiation Interchange Factor 
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C d and Nu stand for drag coefficient and Nusselt number for heat transfer which are 
functions of Reynolds number and relative Mach number. Typical correlations are given in 
Refs. 3.16 and 3.17. Carlson and Hoglund’s correlation (Ref. 3.16) is written as: 

C d = (24/Re) (1 + 0.15 Re 0 687 ) (1 + e-)/ 

[1 + M (3.82 + 1.28 e 1 25Re/M )/Re] 

Nu = (1 + 0.2295 Re 055 )/ [1 + 3.42 M (2 + 0.459 Re 055 )/Re] 

where a = 0.427/M 463 + 3.0/Re 0 88 . A more accurate but more complicated correlation for 
the drag coefficient is provided by Henderson (Ref. 3.17). That is, for Mach > 1, 

C d = 24 [Re + S {4.33 + exp(-0.247 Re/S) (3.65 - 1.53 T„/T) 

/(I + 0.353 T w /T)}] ' 

+ exp(-0.5*M/Re 1/2 )[0. 1M 2 + 0.2M 8 + (4.5 + 0.38a) 

/(I + a)] + 0.6 S [1 - exp(-M/Re)] 

where S = M(y/2) 1 2 is the molecular speed ratio, a = 0.03 Re + 0.48 Re 1/2 . For Mach > 
1.75, 

C d = [0.9 + 0.34/M 2 + 1.86(M/Re) 1/2 {2 + 2/S 2 

+ 1.058 (TJT) m /S - 1/S 4 }] / [1 + 1.86 (M/Re) ,/2 ] 

And, for 1 < Mach < 1.75, 

C d = C d M=1 + (4/3) (M - 1) (C dM=1 . 7S - C dM=1 ) 
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which assumes a linear variation between M = 1 and M = 1.75. 

While the Henderson drag law has been found to give slightly better motor performance 
predictions, the differences in results using the Carlson-Hogland and Henderson method are 
slight. The Henderson method was used in the Cycle 2 RAMP/SPF2 (Ref. 3.5) model and can 
easily be incorporated into FDNS. All FDNS calculations presented in this report used the 
Carlson-Hogland model. The Nusselt number correlation of Drake (Ref. 3.18) which 
corresponds to the Cycle 1.0 methodology was used for all FDNS calculations. It is 
recommended that the heat transfer model of Moylan (Ref. 3.19), which was developed for the 
Cycle 2.0 plume methodology, be incorporated into FDNS. 


Details of the Particle Solution Method . In the present two-phase flow model, an independent 
module was employed for the calculation of particle drag forces and heat flux contributions to 
the gas flow field. Subroutines for locating the particles and integrating their trajectories are 
called for each particle size group. The drag forces and heat fluxes are then saved for every 
grid point. These forces and fluxes are then used to evaluate the particle source terms in the 
gas-phase governing equations. In the present FDNS flow solver, either of two forms of the 
energy equation (i.e. static enthalpy form or total enthalpy form) can be selected. It has been 
found that although either form of energy equation usually gives similar solutions, the static 
enthalpy equation provides better definition of the liquid rocket plume shear layers, as shown 
by extensive solutions made for the SSME. The energy equation presented previously under the 
governing equations section is the total enthalpy form. The static enthalpy option (see Section 
3.2.2) should be used for two-phase flow solutions. 

Particle wall-boundary conditions are treated by using a specified fraction of the colliding 
particles which stick to the wall. Particles which stick result in a decreased particle velocity 
normal to the wall for that particle size fraction. Therefore, for the particle size fraction which 
locally collides with the wall, part of the particles stick and the other part is turned parallel to 
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the wall. Energy exchange is assumed to be due only to the particles which stick. This model 
of particle wall interaction can be improved, but new experimental test data must become 
available in order to do so. 

In the 2-D version of the FDNS flow solver, a fourth-order Runge-Kutta method was 
employed to integrate the particle trajectories. After a thorough test of the integration routine, 
it was found that the explicit scheme sometimes results in divergent particle solutions when the 
source terms become large. Therefore, an implicit integration scheme was employed in the 
present model. For convenience, consider the X-component of the particle equation of motion. 
That is, 

dXp/dt = Up 
dU p /dt = A (U c - Up) 
where A = l/t d 

U c = gas velocity 
Up = particle velocity 
X p = particle location 

In finite difference form the above equations can be written as: 

XpO+D . Xp (n) = (At/2) [U p (n+1) + U p (n> ] 

U p (n+,) - U p (n) = AtA [U c - Up (n+1) ] 
or 

X p (0+,) = Xp (n) + At/2 [U p (0+,) + U p (n) ] 

UpO+n = [u p w + AtA U c ]/(l+AtA) 

These two equations are unconditionally stable despite the magnitude of the source terms. 
To provide better time resolution, a variable time step size is chosen so that a particle would 
take at least 4 time steps to go across a grid cell. 
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3.2.2 History and Validation of the FDNSEL Navier Stokes Two-Phase Code 

The two-phase flow capability was added under a previous NASA study (Ref. 3.7) to 
support the development of a solid rocket motor plume impingement model for predicting launch 
stand environments. Checkout cases for the previous study focused on modeling the flowfield 
of a 20% scale model of the Space Shuttle solid rocket motor combustion chamber/nozzle 
flowfield. Results of the FDNS calculation for this motor were then compared with a RAMP2 
nozzle solution for the same case. While the results of these comparisons were for the most part 
qualitatively acceptable, quantitatively there were enough differences in the results that the 
application of an FDNS flowfield for radiation predictions was not recommended. One of the 
reasons that the results were not absolutely comparable was the geometry which is used in the 
combustion chamber. The combustion chamber for the FDNS calculation was simulated as 
shown in Fig. 3.1 which corresponds to a simulation of the grain geometry late in the bum of 
the motor. The RAMP2 transonic module assumes an infinite sink at an inlet angle 
corresponding to the inlet angle to the throat. The FDNS code was also run for turbulent 
reacting flow while RAMP2 was run using laminar equilibrium chemistry. For these reasons 
it was not possible to absolutely check out the FDNS solution during the previous study. 

At the onset of this study, validation of the FDNS code continued using the same check 
case as was used previously with little improvement in the comparative results. However, 
instead of the radical geometry used in the previous comparison, a more regular geometry as is 
shown in Fig. 3.2 was used. This corresponded to an early bum time. All other variables that 
dictate the solution were identical, i.e. : 

• Frozen chemistry 

• Prandtl number = .7 

• Viscosity and viscosity exponent (.6 laminar) 

• Particle-gas heat transfer model (Drake) 

• Drag law - Carlson-Hogland 
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Figure 3.1 FDNSEL MNAS A08/ ASRM48-5 Nozzle Flowfield Grid (179 x 81) 
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Fig. 3.2 FDNSEL MNASA08 ASRM Contoured Nozzle Grid (201 x 41) Simulating Early Bum Time 
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Results of RAMP2 and FDNS calculations are shown in Figs. 3.3 and 3.4. These figures 
show a comparison of Mach number and temperature distributions at the exit plane of the motor. 
The results are significantly different. It is apparent that the majority of the particles are 
contained in a smaller area of the exit for the FDNS solution. This is evident by observing the 
peak in temperature at .9 ft radius while the RAMP solution peaks at 1.05 ft. This effect could 
be attributed to the difference in the combustion chamber geometry used for FDNS versus that 
used by the RAMP transonic module. Also, notice the spike in temperature that FDNS predicts 
near the axis. This was traced to predicted particle number densities on and near the axis. The 
FDNS calculation used a single particle trajectory at each grid point to perform the Lagrangian 
tracking. This, compounded with the tight grid near the axis, led to numerical problems with 
the code that resulted in a poor distribution of particle number densities. This could have been 
corrected by using more trajectories in each cell and changing the grid; but even then the results 
of the RAMP and FDNS calculations would be different enough that any conclusions about the 
accuracy of FDNS would not be possible. At this stage of the validation, it was decided to 
eliminate combustion chamber geometry effects and concentrate on validating the equations 
which are solved by FDNS for two-phase flow. 

A 15 degree source flow case was set up for RAMP and FDNS. Identical start lines 
were input to both codes consisting of the following conditions: 

• Mach number - 2.0 

• Gas temperature - 6000 

• Gas velocity - 6500 ft/sec 

• Molecular weight - 20 

• Pressure - 500 psi 

• Particle size - 4 micron radius 

• Particle/gas flow rate ratio - .5 

• Particle temperature - 6500 °R 

• Particle velocity - 6000 ft/sec 
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Fig. 3.4 Temperature Profile at the MNASA ASRM Contoured Nozzle Exit 
(Slip Wall, Two-phase, Frozen Chemistry) 
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Calculations using FDNS resulted in the flow going subsonic downstream of the start line 
while RAMP had a slight decrease in Mach number near the start line followed by a gradual 
acceleration of the flow. 

First, it was thought that the initial guess for the flowfield was the problem. To check 
this out the initial flowfield was set to exactly what RAMP calculated. The flow still went 
subsonic. Next, the calculated drag and heat transfer terms were compared to those calculated 
by RAMP. They were found to be the same. By examining the trend in the results, evidence 
pointed toward the gas energy equation since far too much energy was transferred to the gas 
which caused the flow to heat up and decelerate. Upon looking at the terms in the static form 
of the energy equation, it was found that the sign and magnitude of the work loss portion of 
equation was incorrect. Instead of multiplying the difference in gas and particle velocity by the 
drag force, the code was multiplying the absolute particle velocity by the drag force (which is 
the total energy form of the equation). The energy equation was modified and the calculation 
rerun. Figures 3.5 and 3.6 show a comparison of RAMP and FDNS for pressure and 
temperature distributions along any given gas streamline. The results are almost exactly the 
same. To further verify the energy equation, several axisymmetric constant area duct flow cases 
were set up for RAMP and FDNS. 

Constant area, duct flow, two-phase cases eliminate particle trajectory effects since the 
particle streamlines remain straight and particle number density is only affected by the change 
in particle velocity. Several cases were run making various assumptions on particle and gas 
temperature and velocity lags. The results of three of these cases are shown in Figs. 3.7 thru 
3.12. All cases assumed a particle/gas flow rate ratio of 0.5 and a static pressure of 500 psi. 
Case 1 assumed that the particle and gas velocities were 6500 ft/sec, the particle temperature 
was 5500°R and the gas temperature was 5000°R. Figure 3.7 presents a comparison of the gas 
and particle temperature distributions down the duct for FDNS and RAMP. The results are 
almost identical. Differences in RAMP and FDNS at the beginning of the duct are due to 
differences in step size and the fact that FDNS uses gas properties corresponding to the flow 
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properties at the beginning of each particle trajectory step while RAMP uses average gas 
properties. Figure 3.8 presents velocity distributions for Case 1. Case 2 is the same as Case 
1 except the particle velocity is 5500 ft/sec and the gas velocity if 6500 ft/sec. Again, the 
temperature (Fig. 3.9) and velocity (Fig. 3.10) distributions are almost identical. Case 3 is the 
same as Case 1 except the particle velocity is 7500 ft/sec while the gas velocity is 6500 ft/sec. 
Figure 3.11 presents a comparison of the temperature distribution and Fig. 3.12 presents the 
axial distribution of particle and gas velocity. Additional cases making various other 
assumptions on gas and particle temperature and velocity lags were calculated with similar 
results. These results confirm that the momentum and energy transfer between the particles and 
gas are now properly described by the FDNS governing equations. 

It now appears that FDNS is solving the proper set of equations. In order to determine 
the effect and the difference in treatment of the transonic region by FDNS and RAMP, the 
nozzle case was rerun using the corrected version of FDNS for the geometry shown in Fig. 3.2 
as well as a new geometry corresponding to a later bum time which is shown in Fig. 3.13. 

Results of the two FDNS and RAMP calculations are presented in Figs. 3. 14 thru 3.17. 
Figures 3.14 and 3.15 present centerline distributions of Mach number and temperature, 
respectively. The centerline Mach number distributions presented in Fig. 3. 14 show that FDNS 
allows the flow to accelerate more along the centerline of the motor than does RAMP. Inlet 
geometry effects predicted by RAMP shows that for the case where the grain has burned back 
(which results in a steeper effective inlet angle) the flow does not accelerate as much as the 
initial bum case. FDNS results for the two cases show an opposite trend. Centerline 
temperature distributions shown in Fig. 3.15 indicate similar trends with RAMP having higher 
centerline temperatures than FDNS, as well as opposite trends with bumback geometry changes. 
Figures 3.16 and 3. 17 present exit plane Mach number and temperature distributions for the two 
FDNS and RAMP cases. In Fig. 3.16, RAMP predicts less acceleration of the flow except near 
the outer portion of the flow where the particle limiting streamlines are located. The overall 
trends of flow acceleration of the two RAMP cases versus the FDNS cases are again reversed 
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Fig. 3.15 Temperature Distribution Along the Centerline of a Contoured Nozzle 
(Two-Phase, Frozen Chemistry, Laminar) 
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Fig. 3.16 Mach Number Profile at the Exit of a Contoured Nozzle 
(Two-Phase, Frozen Chemistry, Laminar) 
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as was previously shown for the centerline results. Figure 3.17 shows similar trends to the 
Mach number distributions with RAMP2 predicting overall higher temperature. The exit plane 
temperature distributions indicate that the steeper inlet angle results in containing a majority of 
the particles in a smaller area of the exit. This is apparent by observing the location in the exit 
where the temperature starts to rapidly drop off. For the RAMP cases, this occurs at approxi- 
mately 1.05 feet for the bumback case and 0.95 feet for the initial bum case. FDNS shows this 
occurring at 0.9 feet for the bum back case and 0.85 feet for the initial bum case. The 
implications of the results presented in Figs. 3.14 thru 3.17 are that there should be significant 
differences in the particle density distributions between the RAMP and FDNS predictions. 

Two possible explanations for the observed differences in the RAMP2 and FDNS results 
shown in Figs. 3.14 thru 3.15, that would influence the particle density distributions, are the 
Legrangian tracking method used by FDNS and combustion chamber geometry differences. 

At the present time, FDNS assumes that the mass flux of the particles is constant along 
the initial data surface, although particulates need not be present at all points on this surface. 
The user may also specify how many particle trajectories may be initiated for each particle size 
at each cell on the initial data surface. The way the Lagrangian tracking method works is to use 
the trajectory information to effectively determine how many trajectories go through each cell 
and then allocate the mass associated with these trajectories to the particle terms in the forcing 
functions of the gas equations at each of the points that define any particular cell. Typically, 
FDNS is run with one trajectory for each cell. For uniformly expanding cases one particle 
trajectory provides enough accuracy to produce good results which is demonstrated by the source 
flow and duct cases presented earlier in the report. However, for the nozzle case where the flow 
is contracting and expanding, one trajectory may not be adequate. Figure 3.18 presents exit 
plane number density distributions for the smallest and largest particle size at the nozzle exit 
plane for the initial bum geometry nozzle case. Contained on this figure are results for RAMP 
and 3 FDNS cases. The FDNS results are for cases where 1, 5 and 20 particle trajectories were 
initiated in each cell. It is apparent from this figure that the particle number density becomes 
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Fig. 3.18 MNASA ASRM Contoured Nozzle Exit Plane Number Density Distribution Using FDNS and RAMP2 
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better behaved as the number of trajectories increase, especially for particle 1 which is the 
smallest particle. This is to be expected since the smaller particles are more easily influenced 
by the gas. Also, of note in Fig. 3.18 is the large dip in number density near the axis for the 
smallest particle size. This is artificial, and is due to two factors. First, the FDNS code has 
a reflecting boundary condition at the axis of symmetry and walls. Thus, any particles that 
intersect the boundary are reflected, and in theory it is possible that FDNS could predict no 
particles at the centerline. Secondly, if the grid distribution is very fine near the axis, FDNS 
could also predict no particles at the axis. In reality this is not the case. The Lagrangian method 
needs to be improved near boundaries, perhaps with an extrapolation method. The dip in larger 
particle number density near the limiting streamline is caused by particles intersecting the nozzle 
wall and being reflected. More work needs to be done on FDNS in the treatment of particles 
near the boundaries. In spite of the improved number density distributions that resulted from 
using 20 particle trajectories per grid cell, there was little effect on the temperature of the 
flowfield. Two further calculations were made to help with the interpretation of the number 
density results. The number densities near the axis are approximately 20% below those 
predicted by RAMP. To verify that the gas results are consistent with the particle number 
density, RAMP2 was run for the same case but with the particle/gas flow rate ratio reduced by 
20%. The predicted temperature at the exit plane was very close to that predicted by FDNS 
(~3600°R at the exit plane centerline). This further confirms that FDNS is properly handling 
the particle-gas interaction. As a final test of the calculation of particle number density, an 
additional check was made on the number density using another method. 

As part of the Cycle 2 plume methodology (Ref. 3.5), a particle trajectory tracing code 
was developed. This code traces trajectories through a known gas flowfield and can be used to 
calculate number densities if one knows what the particle number density is at the point where 
the particle trajectory is initiated. The FDNS flowfield was mapped in a format that could be 
used by the trajectory program. Particle trajectories were initiated in this flowfield and tracked 
through the mapped flowfield. The predicted trajectories using the trajectory code were almost 
exactly what was predicted by FDNS. Using a streamtube/mass flow balance of the trajectory 
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code results gave number densities near the axis that were consistent with FDNS predictions. 
Therefore, by the process of elimination, the observed differences in RAMP2 and FDNS for the 
combustion chamber/nozzle case is most likely the result of geometric differences in combustion 
chamber treatment of the two codes. 

The transonic module in the RAMP2 code assumes sink flow at the entrance angle to the 
throat (30° for the early bum case). FDNS was set up for the grid shown in Fig. 3.2 which has 
the flow entering a gentle contracting region followed by the 30° entrance to the throat. Both 
RAMP2 and FDNS assume that the particle mass flux is uniform across the inlet so that the 
initial values used by the two codes won’t be responsible for the observed difference in the 
results. The only absolute confirmation that the differences in geometry is the reason for the 
difference in the results would be to rerun FDNS with the same source-like geometry. 
However, this was not done due to the stage in the study effort that deficiencies in the energy 
equation was discovered in FDNS. However, qualitatively the differences in the geometry can 
explain the differences in the results. One would expect the RAMP geometry to direct the 
particles toward the axis of the nozzle, since all particle trajectories and gas streamlines are 
focused toward the axis. On the other hand, in the entrance to the converging section of the 
chamber for the FDNS geometry, all particle trajectories (as well as gas flow) is parallel to the 
axis. One would expect that RAMP2 would predict higher number densities in the vicinity of 
the axis under these conditions. 

RAMP2/SPF2 flowfields have been used extensively in the prediction of radiation 
environments. These calculations have shown excellent correlation with measured date. The 
only attempt at validating the FDNS methodology relative to radiation from the flowfield has 
been the study done under this contract. Until that time when more comparisons can be made 
to absolutely show the effects of combustion chamber nozzle geometry on radiation loads, FDNS 
modeling of the combustion chamber for radiation predictions should utilize the same sink flow 
geometry that RAMP2 uses. It is possible that future studies using FDNS modeling of the 
combustion chamber might explain some observed inconsistencies in radiation distributions 
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radially across the plume. 

3.2.3 FDNS Input Instructions and Sample Case For Two-Phase Nozzle Analysis 

The use and preparation of FDNS input files and common block lengths is described in 
detail in reference 3.11. However, some additional comments can be made for running two- 
phase flow cases. To initiate a two-phase flow case, two parameters must be set in the fdnsOl 
include file. UKPMX must be set to IIQMAX to invoke the two-phase Lagrangian tracking. 
The number of particle trajectories to track for each cell is set using the parameter NPMAX. 
For converging-diverging nozzle flows NPMAX should be set to at least 10. FDNS has two 
options for starting or restarting the calculations. If the start option (IDATA=2) is used, then 
the user must input flow and initial flowfield files using the include file fexmpOl. The other 
option for starting FDNS is initiated by setting IDATA=1. In this case, the grid and flowfield 
files must be input to the code via the fort. 12 (grid) and fort. 13 (flowfield) files. The format 
of these files was previously specified in Tables 2.4 and 2.5 of Section 2.2.3. It is usually 
easier to write a code that will generate these files than it is using the IDATA=2 option. 


Flowfield initialization can be very important in obtaining a converged solution for two- 
phase flow cases. For combustion chamber/nozzle cases the initial guess is not as important 
as for fixed upstream (supersonic) boundary cases. If the initial guess for a supersonic case is 
unrealistic, it is possible that a real solution will not be obtained. For combustion 
chamber/nozzle cases, a solution will be obtained, but if the initial guess is poor, excessive 
computer resources will be required for a solution. For these reasons, it is important that the 
initial flowfield (fort. 13) file be as realistic as possible. 

The other user supplied files is the fort. 11 file that was previously described in Table 
2.7. This file provides the overall control variables, boundary conditions, and reference 
properties for the FDNS solution. 
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The sample case that is described in this section is for the geometry previously shown 
in Fig. 3.2. The geometry corresponds to a simmulation of the MNASA48 contoured ASRM 
geometry at a chamber pressure of 590 psia. The fort. 12 and fort. 13 files for this case are 
contained on MS-DOS disks RAD3, RAD4 and RAD5. The fort. 11 file is contained on disk 
RAD3. The source code for the version of FDNS that was used to run this case as well as the 
other cases presented in the report is found on disk RAD6. 

This test case is an axisymmetric nozzle flowfield simulation with frozen gas phase 
chemistry at t=5 sec. motor bum time. The FDNS input required for this case is shown in 
Table 3.1. This table represents the complete FDNS input file for a 201 axial by 41 radial grid. 
The gaseous specie thermodynamic data (in JANNAF/CEC standard format) is for 12 gas 
species (NGAS = 12). The particle input for 5 A1 2 0 3 particle size classes (IDPTCL) of mass 
diameter (DDPTCL) 2.98 5.16, 7.04, 8.69, and 11.7 fim. Each of the particles has a mass 
density of 188 lbm/ft 3 and a total enthalpy of .642149E+08 ff/sec 2 . This enthalpy corresponds 
to a temperature of 6321.6 °R which is the temperature of the gas at the inlet plane. This 
enthalpy is calculated assuming a specific heat of liquid A1 2 0 3 of .34 BTU/lbm/°R, a solid 
specific heat of .32 BTU/lbm/°R, a melt temperature of 4172.4°R, and a heat of fusion of 499.74 
BTU/lbm. The particles were assumed to be in thermal and dynamic equilibrium with the gas 
(UUPTL=U part /U g „= 1.0). The particle mass flow for each particle class is calculated based on 
the particle to gas flowrate predicted by the RAMP2 code and the percentages for each particle 
class. RAMP2 calculates a gas flowrate of 194.66 and a particle flowrate of 97.666 lbm/sec. 
The distribution of flowrates amongst the particle sizes was selected based on the Cycle 2.0 
methodology described in Ref. 3.5. Input to FDNS for each group is the particulate mass flow 
for that size group divided by 2 t. The particles are assumed to be uniformly distributed from 
the nozzle axis to the nozzle wall thru MPTCL=1 (axis) and MPTCL2=41 (wall). 

In the event that the user wants to consider reacting chemistry, NREACT can be set to 
18 and the reaction set shown in Table 3.2 can be added after the reaction header record. 
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Table 3.1 Listing of fort. 11 Input File For Two-Phase FDNS Sample Case 


IDIMj, (2-D axisy. nozzle test case, multi — phase flow) 


IZON,IZFACE, 

IBND, 

ID, 

ISNGL, 




1, 

0, 

4, 

0, 

0, 




IZT, 

JZT 5 

KZT , 

CBGX , 

CBGY , 

CBGZ , 

CBVX , 

CBVY, 

201 , 

41, 

1, 

0.0, 

0.0, 

0. , 

0., 

0., 

NNBC, 

IZB1 , 

IZF1, 

IJZ1, 

I JZ2 , 

JKZ1 , 

JKZ2 , 

(2 LINi 


IZB2 , 

IZF2 , 

IJZ1, 

I JZ2 , 

JKZ1 , 

JKZ2, 


IBCZON, 

IDBC , 

ITYBC , 

IJBB, 

UBS, 

I JBT , 

I KBS , 

I KBT , 

1, 

A- H 

1, 

1, 

1, 

41, 

1, 

1, 

1, 

4, 

3 , 

1 , 

1, 

201, 

1, 

1, 

1, 

3 , 

3 , 

41, 

1, 

201, 

1, 

1, 

1, 

1, 

o 

201 , 

1, 

41, 

1, 

1, 

IWBZON, 

r. 

H 

u 

2, Ml, 

M2, 

Nl, N 

2, IWTM , 

HQDOX , 

IWALL 

ISNZOM, 

I SNBC , 

ISNAX, 

ISNBS, 

ISMBT , 




I DAT A, 

I GEO , 

6ITT, 

ITPNT, 

I COUP, 

NLIMT, 

I AX, 

ICYC , 

1, 

41 , 

201, 

200, 

3 , 

1, 


0 , 

-5.000E- 

01DTT , 

IREC, 

REC, 

THETA, 

BETAP, 

IEXX , 

PRAT, 

5.000E-01 , 


1.00, 

1.0, 

1.0, 

1, 

-1.0, 

I PC, 

JPC , 

IPEX, 

JPEX , 

IMN, 

JMN, 



202, 

1, 

4221, 

1, 

26o , 

1, 



VISC ( 1/RE ) , IG , 

ITURB , 

AMC, 

GAMA, 

CBE , 

CBH , 

EREXT, 

9.09960E 

-07, 1, 

O 

Jl 

0.089, 

1.265, 

0.0, 

0.0,1 

. E-08 , 

ISWU, 

ISWP, 

ISWK, 

I SKEW , 





1, 

1, 

1, 

0 ? 





INSO(IEQ) : (VISCOSITY 

■ 4.422 

8260E- 

07 SLUGS/FT-SEC) 

U, V, W 

,TN., I)K , 

DE, 7, 

8 , 9, VS 

,FM,SP 

!« 



1, 1, 0 

, 1 , 0, 

0, 0, 

0, 0, 0 

,1,1 

j, 



NGAS , 

NREA18 

, IUNI1 

, DREF(SLG), 

UREF(FZS) 

, TREF(R), XF 


CBVZ , 


DEN NX , VISWX, 


H20 

0.26340654E+01 0 . 31 121899E— 01 


(FT) , 

. .00000E-0 , 

300-000 5000.000* 18.01520 

-0 - 902784 5 IE— 06 0 . 12673054E-09-0 - 69164734E-14 
-0 . 298762 58E+0 5 0 . 70823874E+01 0. 41675563E+01-0 . 18106S68E-02 0 . 594 50877E-05 
-0 . 48670872E-08 0 . 1 5234 144E- 11-0. 30289 547E+05-0 . 73087996E+00 

300.000 5000.000 31.99880 

19820646E-06 0 - 33749007E-10-0 . 23907374E-14 
37837 136E+0 1—0 . 30233634E-02 0 . 99492754E-05 
10638 107E+04 0 . 3641634 5E+01 

300.000 5000.000 2.01580 


D2 

0.36122139E+01 
— 0 .11 9731 51E+04 0 . 36703308E+01 
—0 .981 8910 IE— 08 0 . 33031826E— 1 1 
H2 


0 . 74853166E-03-0 . 

0 . 

- 0 . 


0.30558124E+01 0 . 59740403E-03-0 . 16747471E-08-0 . 21247544E-10 0 . 25195486E-14 


-0.86168475E+03-0. 17207073E+01 0.29431 


28E+01 0 . 3481 5508E-82-0 . 77713821E-05 


0 .74997493E-08-0 . 25203379E-1 1-0 . 976954 10E+03-0 . 13186136E+01 
0 300.000 5000.000 15-99940 

0.25342960EH-91— 0 . 12473170E— 04-0 . 12562724E -07 0.69029860E-11-0.63797098E-15 
0. 29231 107E+05 0 . 4962B592E+01 0.30309401E+01-0.22525853E-02 0 . 39824 540E-05 
— 0 . 3260492 IE— 08 0 . 10152035E-1 1 0 . 29136525E+05 0 . 26099341E+01 


H 


300.000 5000.000 


1.00790 


0.25000000E-+-01 0-00000000E+00 0.00000000E-+-00 0.00000000E+00 0 . 00000000E+00 


0 . 2 547439 1E-+-0 5—0 . 4 598984 1E+00 0. 
0.00000000E+-00 0.00000000E+00 0. 


?5000000E+01 0 - 00000000E+00 0.0Q000000E+00 
2547439 1E-+-0 5—0 . 4 598984 1E+00 

300.000 5000.000 17.00730 

0 . 22048808E— 06 0 .20191288E-10-0.39409830E-1 5 
0 . 38737299E+01— 0 . 13393773E-02 0 . 16348351E-05 
0 . 35802349E+04 0 . 34202406E+00 

300.000 5000.000 28.01040 

0.29840696E+01 0 . 14891390E-02-0 . 57899683E-06 0. 10364577E-09-0 . 69353 550E-14 
— 0 . 1424 5228E+05 0 . 634791 56E+01 0 - 371 00928E+01 -0 . 16190964E-02 0 . 36923593E-05 
—0.203 19675E— 08 0 ■ 23953344E-12-0 . 143 5631 0E+05 0 . 29555352E+01 

300.000 5000.000 44.00980 

0.30981719E-02— 0. 12392571E-05 0 . 22741325E-09-0 . 1 5525955E-13 
-0 . 48961441 E+0 5—0 . 98635983E+00 0 . 24007797E+01 0 . 87350961 E-02-0 . 66070879E-05 
0 . 20021362E-08 0 . 63274039E-1 5-0 . 48377527E+05 0 . 96951456E+01 


OH 

0. 288978 15E-+-01 0 . 10005879E-02- 
0 . 388 5704 1E+04 0 . 5556642 5E+01 
— 0 . 52133636E— 09 0 . 41826975E-13 
CO 


C02 

0.44608040E+01 
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Table 3.1 Listing of fort. 11 Input File For Two-Phase FDNS Sample Case (Cont.) 


CL??? 300.000 5000.000 35. 

0.29537796E+01-0.40792712E-03 0. 1 5288342E-06-0 .2638434 5E-10 0.17206581E- 
0- 13695677E+05 0. 30667325E+01 0 . 20774281E+01 0 . 29487 169E-02-0 . 4391 9732E- 
0 . 24499776E-08— 0 . 41007685E— 12 0 . 13871928E+05 0 . 73136343E+01 
CL2?? 300.000 5000.000 70. 

0 - 4307781 4E-+-01 0 .31 182816E-03-0 . 16310807E-06 0 .4451 1913E-10-0 . 43057753E- 
-0.13458251E+04 0 . 206666S4E+01 0 . 31316886E+01 0.48997877E-02-0. 6941 1463E- 
0 . 4478 564 IE— 08-0 . 10621859E-1 1-0 . 10979696E+04 0 - 77833424E+01 
HCL?? 300.000 5000.000 36. 

0 . 27665884E-+-01 0 . 14381883E-02-0 . 46993000E-06 0 . 73499408E-10-0 . 4373 1 106E- 
-0. 11917468E+05 0 . 64583540E-M31 0 . 3524Q171E+01 0 . 29984862E-04-0 . 86221891E- 
0 . 2097972 IE— 08-0 . 98658191E— 12— 0 .121 50509E+05 0 . 2395771oE+01 
N2 300.000 5000.000 28. 

0. 28532898E+01 0 . 16022128E-02-0 . 62936891E-06 0. 11441022E-09-0.78057466E- 
— 0 . 89008093E+03 0 . 63964896E+01 0 . 37044177E+01 -0 . 14218753E-02 0 . 28670393E- 
-0 . 1 202888 5E -08-0 . 13954677E-13-0 . 10640795E+04 0 . 22336285E+01 
IDPTCL.,5*** PARTICLE INPUT CONTROL **** 


45300 

-14 

■05 

90600 

-14 

•05 

46100 

■14 

•06 

01340 

■14 

05 


5, 


0 , 


I F'TZON 5 I DBCPT , 

LPTCL1 * LPTCL2 , 

flPTCLl „ 

MPTCL2 p NPTCL 1 , NPTCU 

ITPTCL, 

DDPTCL , DNPTCL , 

WDMASS , 

UUPTCL , HTPTCL , 

1 , Ip 

1 , 2 , 

1 , 

41, 1, 

5 , 

2-980,188.00, 

0.307, 

1 . 000 , 0 . 642149E-+-08 

1 , 1 , 

Ip 2 , 

li 

41, 1 , 

5, 

5.160,188.00, 

1 - 80S 

1 . 000 , 0 . 642149E+08 

1 , 1 , 

Ip 2 , 

1 , 

41, 1, 

5, 

7.040,188.00, 

4.310, 

1 . 000 , 0 . 642149E+08 

1 , 1 , 

1 j, 2 , 

1 , 

41, 1, 

3, 

8.690,188.00, 

o . 662 , 

1 . 000 , 0 . 642149E+08 

1 , 1 , 

Ip 2 , 

1 , 

41, 1, 

5, 

11- 70 , 188 - 00 , 

5.457, 

1 . 000 , 0 . 642149E+08 


LINES EACH) 


1 , 
3, 
1 , 
3 i 
1 , 
3 , 

3 , 
1 , 


***END****** 
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Table 3.2 Listing of FDNS Finite Rate Chemistry Reaction Set For fort. 11 Input File 


REACTION: H20, 02, H2 

, 0 , H , 

OH, 

CO, 

C02, 

CL. 

CL2, 

HCL.N2, 

1, 

1.7000E13, 0.00, 

24070., 0, 

0. 








e., -i., -i. 

, 0. , 0. , 

n 

0- 5 

0-5 

0- , 

0- 9 

0. 

90- 9 

2, 

2 . 1900E13 , 0.00, 

2590., 0, 

0. 








1., 0., -1. 

, 0., 1., 

“1-5 

0-5 

0. , 

0-9 

0-9 

0. 

90- 9 

3, 

6.0230E12, 0.00, 

550., 0, 

05 








1., 0., 0. 

, 1., 0., 

~ - , 

0. , 

0-9 

0-9 

0-9 

0. 

9 0- 9 

4, 

1.8000E10. -1.00, 

4480., 0, 

05 








0., 0., -1. 

, -1. , 1. , 

1-9 

0- , 

0. , 

0- , 

0-9 

0. 

90- 9 


1.2200E17, 0.91, 

8369., 0, 

09 








0., -1., 0. 

, 1. , -1. , 

1-9 

0- 5 

0-9 

0-9 

0- 9 

0. 

,0- , 


4 . 0000E12 , 0.00, 

4030., 0, 

09 








0., 0., 0. 

, 0. , 1. , 

"1- , 

“1-5 

1-9 

0. , 

0-9 

0. 

,0-9 

7, 

3 . 0000E12 , 0.00, 

25000. , 0, 

0 1 








0., -1., 0. 

, I-, 

0- 9 

“1-5 

1-9 

0- 9 

0-9 

0. 

,0- 9 

8., 

1 -0000E16 , 0.00, 

0 . , 999 , 

05 








0., 0., 0. 

, -1., -1., 

1-5 

0- , 

0- 9 

0. 9 

0- 9 

0. 

,0. , 

9, 

2.5500E18, 1.00, 

59390. ,999, 

0, 








0. , 1 . , 0. 

, -“2. , 0. , 

0. , 

0. , 

B. , 

6. , 

Q. , 

0. 

,0- , 

10, 

5 . 0000E1 5 , 0.00, 

0- ,999, 

05 








0. , 0. , 1 . 

, 0., -2., 

0. , 

0- 5 

0-9 

0. , 

0. , 

0. 

,0. , 

11, 

8.4000E21 , 2.00, 

0- ,999, 

0, 








1., 0., 0. 

, 0., -1., 

~l-9 

0- 5 

0-9 

0. , 

0-9 

0. 

0. , 

12, 

6.0000E13, 0.00, 

0. ,999, 

0. 








0 . , 0 . , 0 . 

, ~1. , 0. , 

0- 9 

“1-5 

1-9 

0- 9 

0-9 

0. 

,0 - , 

13, 

8.4300E13, 0.00, 

2144., 0, 

0 . 








0., 0., -1. 

, 0., 1., 

B. 9 

0-9 

0. , 

“1-9 

0-9 

1 . 

,0. , 

14, 

3 -0100E13 , 0.00, 

8858., 0, 

09 








-1., 0., 0. 

, 0 . , 0 . , 

1-5 

0. , 

0- 9 

“1-9 

0-9 

1 . 

,0.9 

15, 

3 . 6100E12 , 0.00, 

3020., 0, 

05 








0., 0., 0. 

9 -I- , 0-9 

1 - 9 

0-9 

0. , 

1-9 

0. , 

“1 . 

,0. , 

16, 

9 . 0300E13 , 0.00, 

604., 0, 

0. 








0 . , 0 . , 0 . 

9 0- 9 “1 - 9 

0- 9 

0. , 

0. , 

1-9 

“1-9 

1 . 

,0- , 

17, 

3 . 6300E14 , 0.00, 

~906. ,999, 

05 








0 . , 0. , 0 1 

9 0., 0 . , 

0. , 

0- 5 

0. , 


1-9 

0. 

,0. , 

18, 

1.4500E22, 2.00, 

0. ,999. 

09 








0., 0., 0. 

9 0-9 -1-9 

0- 5 

0-9 

0. , 

“1-9 

0. , 

1 . 

,0. , 
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3.2.4 Influence of Combustion Chamber Geometry on Flowfield Radiation Properties 

One of the observations that has been made based on ground firings of solid rocket 
motors is that there is a tendency for measured radiation to increase with bum time, even if the 
chamber pressure is relatively constant. The MNASA series of tests showed 20-30% difference 
in radiation measurements at the same chamber pressure but at different times in the bum. Two 
possible reasons are carbon due to burning insulation or flowfield effects due to changes in the 
combustion chamber geometry as the grain bums back. Sambamuthi (Ref. 3.20) presented a 
good argument that burning insulation could account for increased radiation with time of bum. 
The results that were previously presented for the MNASA-ASRM combustion chamber/nozzle 
cases suggest that differences in the combustion chamber geometry could result in particle 
number density distribution changes within the nozzle, that in turn could result in different 
particle and gas temperature distributions in the plume. Figure 3. 19 presents SIRRM (Ref. 3.21) 
single line-of-sight radiance calculations at the exit plane for two RAMP2 cases and one FDNS 
case. The two RAMP2 cases correspond to the 30 degree inlet (initial bum) and 54 degree inlet 
(bum back near the end of firing) MNASA contoured ASAM nozzle case. The FDNS results 
are for the 30 degree inlet case shown in Fig. 3.2. Examination of the RAMP2 results show 
that the integrated heat flux for the 54 degree inlet case is approximately 20% higher than the 
initial bum back case (30 degree inlet). These results are consistent with the observed 
differences in the measurements early and late in the bum. FDNS results for the initial bum 
back case are approximately 12% lower than the corresponding RAMP2 case. The absolute 
magnitude of the results shown in Fig. 3.19 are not important since the calculations were 
performed using a Cycle 1.0 methodology with frozen chemistry. The important finding is that 
combustion chamber geometry can influence predicted radiation. Further calculations need to 
be performed to validate and quantify combustion chamber effects. 

3.3 Particle Size Distribution 

One of the uncertainties in performing a two-phase flow calculation is the mean particle 
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size and distribution of mass about the mean size. A majority of the mass of 
aluminum/aluminum oxide particulates that form at the propellant surface are comprised of very 
large particles (>100 fim ) that subsequently break up in the transonic regions of the nozzle due 
to the large shear stresses that are present in this region. To investigate the validity of the 
existing particle size correlations, this portion of the overall study effort used the One- 
Dimensional Three-Phase Reacting Flow with Mass Transfer Between Phases code (OD3P) (Ref. 
3.22). OD3P has been applied to the problem of predicting A1 2 0 3 particle size measured during 
the NASA/MSFC 48 in. diameter ASRM/RSRM subscale solid motor MNASA test series. The 
particle size measurements taken during the MNASA test series are described in Ref. 3.23. In 
the following analysis the mass mean averaged particle diameter (D 43 ) for the 
MNASA9(RSRM48-3) test was predicted using a modified version of the OD3P program. 

The MNASA9(RSRM48-3) test was chosen for analysis because in this test the largest 
number of samples were taken and reported (7 samples) during the test series. The operating 
characteristics at t=5 sec for the RSRM48-3 motor are shown in Table 3.3. The RSRM 
propellant contained: 69.7% Ammonium Perchlorate (AP), composed of 70% by mass of 200 
pm and 30% by mass of 20 pm diameter AP particles; 16% Aluminum particles of 30 /xm 
diameter; 14% PBAN Binder; and 0.30% Iron Oxide. The formulation of the RSRM 
propellant is reported Ref. 3.23. 


Table 3.3. MNASA 09/RSRM48-3 
Motor Operating Characteristics @ t = 5 sec 


Propellant Binder Type 

PBAN 

Aluminum Loading (%) 

16. 

Nominal Chamber Pressure (psia) 

680. 

Calculated Gas Flowrate (lb m/sec) 

241.52 

Calculated Particle Flowrate (lbm/sec) 

96.80 

Throat Diameter, Initial (in) 

9.950 

Throat Diameter, Final (in) 

10.399 

Nozzle Exit Diameter, Initial (in) 

24.104 

Nozzle Type/Liner Material 

Conical/CCP 
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The 0D3P code as documented in Ref. 3.22 and obtained from the Air Force, Phillips 
Laboratory (AFSC), in Aug. 1991 has been reviewed and several improvements have been 
suggested by Mark Salita. The improvements to OD3P suggested by Salita are documented in 
Ref. 3.24, however an updated version of the code including the suggested improvements is not 
available from the Air Force. The work by Salita cited in Ref. 3.24 concluded that one of the 
major deficiencies in OD3P is the model for collision/ coalescence efficiency, particularly 
incomplete particle coalescence efficiency as large particles collide among themselves during the 
contraction/expansion process. 

Based on the OD3P results presented in Ref. 3.24 for the full scale RSRM nozzle which 
indicate that large/large particle collisions have a low coalescence efficiency (approximately 3%), 
and assuming 3 particle sets with a combustion chamber mass median particle diameter of 100 
/<m which can be reduced in size by a factor of 2 upon breakup, the following OD3P particle 
size prediction was obtained for the RSRM48-3 motor. The initial combustion chamber particle 
mass mean diameter of 100 /irn follows the full scale RSRM simulation of Salita in Ref. 3.24 
(115 n m), and the work of Netzer in Ref. 3.25 (approximately 130 jtm near the propellant 
surface). 

The code was modified to assume a constant coalescence efficiency of 3 % , and the input 
particle breakup radius ratio was set to 2.0. The OD3P calculation was initiated in the motor 
chamber at an area ratio of 6.62. The resulting axial particle size predictions are shown in Figs. 
3.20 and 3.21. The OD3P prediction for the three particle group diameters in microns versus 
axial distance normalized by the throat radius (x/r*) is shown in Fig. 3.20. The OD3P code 
predicts that each of the three particle groups will break up three times and reach a final particle 
group diameter of 5.6, 10.8, and 9.9 /*m; starting with initial diameters of 53.7, 100.0, and 
186.1 j«m, respectively. The OD3P prediction for particle mass mean diameter (D 43 ) for all 
three particle groups is compared to the industry standard Hermsen correlation (Ref. 3.26) and 
to the arithmetic average of the seven samples collected during the MNASA/RSRM48-3 test 
reported in Ref. 3.23 in Fig. 3.21. The nozzle exit plane D 43 predicted by OD3P (8.96 ^m) 
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agrees well with the average of the test data for this test (8.64 /xm), however both of these are 
slightly higher (24%) than the Hermsen correlation value of 7.12 /xm, but within the accuracy 
of the correlation. The experimental particle collection technique sampled particles in the nozzle 
exhaust plume from approximately 150 to 620 nozzle exit radii downstream of the nozzle exit 
plane. The assumption is that once the particles have solidified near the nozzle exit plane there 
is no further particle size change downstream of the nozzle exit plane. 

To test the assumption that the particle size does not significantly change beyond the 
nozzle exit plane, the previous analysis of the RSRM48-3 motor has been extended downstream 
beyond the nozzle exit plane in the near field plume to approximately 7 nozzle exit radii. In this 
analysis the one-dimensional, pressure defined, constant area streamtube flowfield option of the 
OD3P program was used to determine the exhaust plume gas and particle properties. The 
analysis was initiated at the nozzle exit plane with gas and particle properties as defined by the 
OD3P nozzle solution. The RSRM48-3 nozzle exit plane gas and particle starting conditions for 
the OD3P plume calculation are noted in Table 3.4. The plume centerline gas axial static 
pressure schedule required as input for the OD3P pressure defined option was determined from 
a typical SPF-II code (Ref. 3. 1) plume sea level flowfield calculation for the RSRM48-3 motor. 
The SPF-II plume centerline axial pressure schedule was normalized by the SPF-II predicted 
centerline pressure at the nozzle exit plane and the result ratioed by the exit plane gas pressure 
as predicted by OD3P which is shown in Table 3.4. In other words the SPF-II code was used 
to determine the shape of the axial pressure schedule, and the initial pressure magnitude was 
determined by the OD3P nozzle solution exit plane pressure. 
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Table 3.4. MNASA09/RSRM48-3 
OD3P Predicted Nozzle Exit Plane Gas and Particle Properties 

at = 5 sec 

Gas Pressure (psia) 

18.73 

Gas Temperature (°R) 

3478.0 

Gas Velocity (ft/sec) 

7687.5 

Gas Density (lb m/ft 3 ) 

0.01007 

| Particle Temperature/Phase || 

• Particle Group 1 (°R) 

3917.9 (subcooling) 

• Particle Group 2 (°R) 

4301.8 (liquid) 

• Particle Group 3 (°R) 

4249.5 (liquid) 

|| Particle Velocity jj 

• Particle Group 1 (ft/sec) 

7219.3 

• Particle Group 2 (ft/sec) 

6726.7 

• Particle Group 3 (ft/sec) 

6796.9 

[| Particle Diameter 

• Particle Group 1 (microns) 

5.63 

• Particle Group 2 (microns) 

10.78 

• Particle Group 3 (microns) 

9.96 

Particle Mass Averaged 
Diameter, D 43 (microns) 

8.96 
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The RSRM48-3 plume centerline axial pressure schedule described above input to the 
OD3P code for this case is shown in Fig. 3.22. The key feature in the predicted plume 
centerline axial gas pressure distribution is the reflected shock which is shown in Fig. 3.22 to 
occur at approximately 4.0 nozzle exit radii. This is typical for a sea level conical nozzle 
expansion for a motor of this type. The plume centerline gas temperature predicted by OD3P 
for this case is shown in Fig. 3.23. The predicted plume particle size axial variation for each 
of the three particle size groups is shown in Fig. 3.24. Particle group 1, the smallest particle 
size group, begins to subcool within the nozzle at a nozzle area ratio of 4.0 (the nozzle exit 
plane is at an area ratio of 5.87); and has completed subcooling and begins to solidify at 2.0 
nozzle exit radii, and has completed solidification at 4.0 nozzle exit radii downstream of the 
nozzle exit plane. Particle group 3, the next largest particle size group, begins to subcool at 0.2 
nozzle exit radii and does not complete subcooling by 7.0 exit radii. Particle group 2, the 
largest particle size group, begins to subcool at 0.4 exit radii and does not complete subcooling 
by 7.0 exit radii. 

The particle mass average diameter versus axial position predicted by OD3P for the 
RSRM48-3 plume is shown in Fig. 3.25. In Fig. 3.25 the OD3P prediction is compared with 
the NASA/MSFC measurements (Ref. 3.23) and the Hermsen correlation (Ref. 3.26) at the 
nozzle exit plane. The prediction is approximately 4% higher than the arithmetic average of the 
measurements for this motor. The most significant finding of this analysis is that once the 
particles reach the motor exit plane and have begun the process of subcooling and subsequently 
solidification, there is little further size change as the particles enter the plume at least to 
approximately 7 exit radii downstream. For this case the particle mass averaged diameter 
decreased from 8.96 /xm at the nozzle exit plane to 8.63 /xm at 6.9 exit radii or 3.7%, which 
is negligible. 

The mean particle size (D 43 ) predicted by OD3P for the MNASA case is very close to 
that measured by Sambamurthi and reported in Ref. 3.23. It is recommended that the 
distribution of particles presented in Ref. 3.23 be used for motors of the size of the MNASA test 
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rsrm48-3 - p 1 ume - od3p axis pressure defined - seca 
center! in© gas pressure 


Fig. 3.22 0D3P RSRM48-3 Plume Centerline Gas Pressure 
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rsrm48~3 — plume — od3p axis pressure defined — seca 
centerline gas temperature 


Fig. 3.23 0D3P RSRM48-3 Plume Centerline Gas Temperature 
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rsrm48-3 - plume - od3p axis pressure defined — seca 
particle d I am© te r C m I c ro ns ) vs. axial position 


Fig. 3.24 0D3P RSRM48-3 Plume Particle Group Diameter vs. Axial Position 

(X/Rj) 
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rsrm48~3 - plume - od3p axis pressure defined - seca 
particle mass avg . dia. (microns) vs. axial position 


Fig. 3.25 0D3P RSRM48-3 Plume Particle Mass Average Diameter vs. Axia 
Position (X/Rj) 
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series and smaller. This distribution was used for Cycle 2.0 predictions. Particle size 
measurements were also taken by Sambamurthi for the full scale RSRM motor (Ref. 3.27). 
These measurements are in very close agreement with the predictions made by Salita using 
OD3P (Ref. 3.24). The mean particle size measured and predicted are almost identical to those 
predicted using the Hermsen correlation. The measured and predicted size and mass fraction 
distributions are again almost identical. The measured mass distribution is best described as a 
monomodal log-normal distribution with a standard of deviation of .13. The results of OD3P 
calculations for the MNASA and full scale RSRM motors when compared with measured 
distributions show that OD3P can satisfactorily be used to make a prior prediction of particle 
sizes and distributions in the plume at least for these two classes of motors. 
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4.0 INTERPRETATION OF MNASA THERMAL IMAGE CAMERA DATA 

The 48 in. diameter MNASA solid motor test series provided diverse and detailed 
measurements of the radiation properties of the exhaust plume that included: total radiometer, 
CVF spectrometer and thermal image camera data. The majority of previous studies used the 
radiometer and CVF spectrometer data to support flowfield and radiation model development. 
This section details the results of a study to investigate the thermal image data relative to the 
radiative properties of solid rocket motor exhaust plumes. The analyses of these additional 
thermal image data sets are required to provide insight into the continuing task of identifying the 
source and magnitude of discrepancies in the comparison of predicted solid motor radiance and 
radiant heat flux vs. measurements. In the discussion presented below, thermal image camera 
radiance data are compared to SIRRM-II code predictions made using both the FDNSEL and 
RAMP/SPF2 flowfield methodologies for various MNASA test series RSRM48 and ASRM48 
subscale motor firings. The intent is to first identify the available MNASA thermal image test 
data, and second to compare radiance predictions to the available thermal image data, and lastly 
to critique the flowfield methodology and to identify shortcomings in the predictive 
methodology. 

The results presented in this section were previously reported in a quarterly progress 
report (Ref. 4.1) and reflect flowfield calculations that were at various levels of maturity and 
validation. The RAMP/SPF2 flowfield solutions are referred to as Cycle 1.0, Cycle 1.5 and 
Cycle 2.0 plumes. These plume calculations refer to three levels of flowfield development that 
resulted from studies (Ref. 4.2) to predict the plume induced environments for the Space Shuttle 
Vehicle equipped with the Advanced Solid Rocket Boosters (ASRB’s). These models are 
summarized in Ref. 4.2. 
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Early in the ASRM plume model development, it was found that the industry standard 
model under predicted the radiation data. One potential hypothesis was that small amounts of 
carbon particulates were present that could potentially result in radiation enhancement that might 
correlate the predictions with measurements. Early comparisons made assuming 2.5% carbon 
by weight in the plume resulted in excellent comparisons with the spectrometer data. Some of 
these comparisons have been included in this section for completeness. However, it was found 
during the ASRM plume studies that carbon was most probably not the factor that caused the 
discrepancy in theory and data. Cycle 2.0 model development identified deficiencies in the 
particle gas heat transfer model and the particle size model that would account for observed 
differences between Cycle 1.0 predictions and data. As a result, Cycle 2.0 plume model is 
recommended for use in predicting solid rocket motor plumes. 

The FDNS results which are presented in this section were made prior to correcting the 
deficiencies in FDNS described in Section 3.2. The main problem with the calculation (i.e. the 
work loss term in the gas energy equation) results in FDNS calculations that simulate a reduced 
heat transfer between the gas and particles. FDNS calculations using the Cycle 2 methodology 
should result in plumes that are similar to those presented in this section for FDNS, albeit for 
different reasons. Due to the stage in the study that this deficiency was identified, it was not 
possible to recalculate the FDNS correlations that are presented in this section. It is 
recommended that these calculations be redone using the corrected version of FDNS. 

4.1 MNASA Test Series Thermal Image Camera Data taken by Sverdrup, Inc. 

Thermal Image Camera data was taken by Sverdrup, Inc. of Arnold Air Force Base, 
Tenn. by V. A. Zaccardi, et al (Ref. 4.3) during the MNASA test series. The data of interest 
to this investigation which have been requested from Sverdrup, Inc./AEDC, are identified in 
Table 4.1. In Table 4.1 the MNASA test number and test date, instrument description, 
bandwidth, and time frame are noted. 
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Table 4.1. MNASA Test Series Thermal Image Camera Data Taken by Sverdrup, Inc. 


Time 

Frame 

(sec) 



Instrument Description 

Bandwidth 

0*m) 

Thermovision Infrared Raster Scanning 
Radiometers: 

• AGA 782-3 

• AGA 782-2 

2.11 to 2.46 
3.14 to 4.08 

AGA 680 Thermovision Infrared Raster 
Scanning Radiometer 

3.41 to 4.00 

Mitsubishi IR-5120AII 
Thermal Image Camera 

2.23 to 2.32 

AGA 680 Thermovision Scanner 

3.41 to 4.00 

AGA 680 Thermovision Scanner 

3.41 to 4.00 

Mitsubishi IR-5120AII 
Thermal Image Camera 

2.23 to 2.32 


5 sec 

9 
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4.1.1 MNASA04/RSRM48-2 Thermovision Thermal Image Camera Data 

The AGA 680 Thermovision Infrared Raster scanning radiometer was used to acquire 
plume thermal images in the 3.41 to 4.00 /xm bandwidth. These data at an unknown time slice 
(but assumed to be between 1 and 27 sec) are shown in Fig. 4.1. The SIRRM-II predicted 
radiance contour plot in the 3.41 to 4.00 /xm bandwidth for the RSRM48-2 test using an SPF2 
flowfield plus 2.5% carbon is shown in Fig. 4.2. The resolution of the thermal image data 
shown in Fig. 4. 1 is not sufficient to make a conclusive comparison with the radiance prediction 
except that the predicted radiance is qualitatively in the ball park, and the maximum predicted 
radiance near the plume centerline at the nozzle exit plane of 0.721 watts/cm 2 /sr is 
approximately 42% of the maximum measured value of 1.2497 watts/cm 2 / sr. 

The Mitsubishi IR-5120AII thermal image camera was used to acquire plume thermal 
images in the 2.23 to 2.32 /xm bandwidth. These data, also at an unknown time slice, are shown 
in Fig. 4.3. The SIRRM-II predicted radiance contour plot in the 2.23 to 2.32 /xm bandwidth 
is shown in Fig. 4.4 for comparison with the measurement. In the comparison of this prediction 
to the measurement, the maximum predicted radiance of 0.350 watts/cm 2 / sr near the plume 
centerline at the nozzle exit plane compares within 9.7% of the measured maximum value of 
0.3875 watts/cm 2 / sr at the same location. If the 2.5% carbon that was added to the SPF2 cycle 

I flowfield is removed, as would be the case in a standard SPF2 flowfield, the resulting SIRRM- 

II radiance contour prediction is shown in Fig. 4.5. In this case the maximum predicted 
radiance of 0.184 watts/cm 2 /sr is a factor of 2.1 below the measured maximum. The 
comparison of the SIRRM-II predicted spectral radiance to the CVF spectrometer measurement 
for wavelengths from 0.7 to 5.7 /xm at a plume centerline location of 0.686 m (or 2.24 nozzle 
exit radii) downstream of the nozzle exit plane for SPF2 cycle 1 flowfields with and without 
2.5% carbon has been shown in Ref. 4.4. This figure, Fig. 4.6, is included here to demonstrate 
that if the flowfield used in the SIRRM-II plume radiance simulation produces a reasonable 
match with the spectrometer data for at least one axial position then the predicted magnitude of 
the maximum radiance will compare reasonably well with the maximum measured radiance as 
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Fig. 4.2 MNASA/04/RSRM48-2 SIRRM-II Predicted Radiance Contours, 3.41 to 4.00 /*m, SPF2 Flowfield plus 2.5% Carbon, 
at t=5 sec 
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Fig. 4.3 MNASA4/RSRM48-2 Mitsubishi IR-5120AII Thermal Image Camera Data, 2.33 to 2.32 urn 
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MNAS A02/RSRM48- 1 SIRRM-II Predicted Plume Radiance and CVC Spectrometer Data, SPF2 Cycle 1 Flowfield, 
x/l — 2.24, at t=5 sec 
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shown in the comparison of Figs. 4.3 and 4.4. 

4.1.2 MNASA03/ASRM48-1 Thermovision Thermal Image Camera Data 

The AGA 782-3 and 782-2 Thermovision infrared raster scanning radiometers were used 
to acquire plume thermal images in the 2.11 to 2.46 and 3.14 to 4.08 ^tm bandwidths, 
respectively. The images from the AGA 782-2 scanning radiometer were saturated at motor 
ignition and therefore not usable for this analysis. The isoradiance data acquired by the AGA 
782-3 scanning radiometer at t = 5.02 sec is shown in Fig. 4.7. In Fig. 4.7, which is a black 
and white representation of the radiance field, it is difficult to determine contour levels and the 
location of the maximum radiance level. However, the maximum measured radiance level 
appears to occur near the plume centerline at the nozzle exit plane and then the plume centerline 
radiance decreases roughly linearly downstream of the nozzle exit plane. 

The SIRRM-II predicted radiance contour plot in the 2.11 to 2.46 /x m bandwidth for the 
ASRM48-1 test at t = 5 sec is shown in Fig. 4.8. In this analysis of the ASRM48-1, the 
flowfield input to the SIRRM-II radiation code was generated using the FDNSEL two-phase 
Navier Stokes flow solver. The FDNSEL flowfield for the ASRM48 motors has been described 
previously in Ref. 4.5, and is used here in preference to an SPF2 flowfield because the FDNSEL 
flowfield best matches the CVF spectrometer data for the ASRM48 test series. Comparison of 
the measured and predicted radiance in Figs. 4.7 and 4.8, respectively, reveals that the predicted 
maximum radiance level of 2.55 watts/cm 2 /sr is 12.1% higher than the measured maximum 
value of 2.275 watts/cm 2 /sr on the plume centerline at the nozzle exit plane; but is located on 
the plume centerline between 1.60 and 2.00 m (5.1 nozzle exit radii) downstream of the nozzle 
exit plane. This characteristic of the ASRM48 motors will be demonstrated in the comparison 
of predicted to measured radiance for other ASRM48 motors presented later in this section. 
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4.2 MNASA Test Series Thermal Image Camera Data Taken by NASA/MSFC 

Thermal image camera data was taken by Don Bryan(ED64) of NASA/MSFC using the 
Inframetrics 600 thermal image camera for tests MNASA04 through MNASA12, excluding 
MNASA09 (Ref. 7). These data for MNASA04 and MNASA06 and a bandwidth of 8 to 12 ptm 
for a single image which is the average of frames from 5.0 to 5.5 sec in time are presented 
herein. 

4.2. 1 MNAS A04/RSRM48-2 Inframetrics Thermal Image Camera Data 

The Inframetrics 600 thermal image camera was used to obtain plume thermal images 
in the 8 to 12 /xm bandwidth for the RSRM48-2 test. These data were taken by Don Bryan 
(ED64) of NASA/MSFC (Ref. 4.6) and provided to us in digital format. Good resolution of the 
visual images has been obtained by averaging several time slices of data together to produce a 
single composite time slice. In this case however, a single frame at the 5.0 sec time slice is 
used. Plume temperature (°F) and radiance (watts/cm 2 /sr) contours from the Inframetrics 600 
camera at the 5 sec single time frame for the RSRM48-2 test are shown in Figs. 4.9 and 4.10, 
respectively. The temperature contours shown in Fig. 4.9 assume an emissivity of 1.0 which 
is a questionable assumption for the assessment of the temperature field of two phase solid motor 
exhaust plumes. The radiance contours shown in Fig. 4.10 were reduced from the actual 
radiance measured by the Inframetrics 600 camera. The instrument was calibrated to a 
maximum temperature of 1345°C (2453°F) and a linear relationship is assumed above this 
temperature. 

The SIRRM-II predicted plume radiance contour plot for the RSRM48-2 using an SPF2 
flowfield plus 2.5% carbon is shown in Fig. 4.11. Comparing the predicted radiance contours 
in Fig. 4.11 with the thermal image measurements in Fig. 4.10 reveals that the predicted 
maximum plume centerline radiance at the nozzle exit plane of 0.557 watts/cm 2 /sr is 7.1% 
higher than the measurement maximum of 0.52 watts/cm 2 /sr at the same location. 
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Fig. 4.8 MNAS A03/ ASRM48- 1 SIRRM-II Predicted Radiance Contours, 2.11 to 2.46 /on, FDNSEL Flowfield, at t= 
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Fig. 4.9 MNASA04/RSRM48-2 Inframetrics 600 Thermal Image Camera Temperature Contours (°F), 8 to 12 pm, 
Emissivity = 1.0, single frame at t=5.0 sec 
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Fig. 4. 10 NNASA04/RSRM48-2 Inframetrics 600 Thermal Image Camera Radiance 
Contours (watts/cm 2 /sr), 8 to 12 /zm, single frame at t=5.0 sec 
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MNAS A04/RSRM48-2 , SIRRM-II Predicted Radiance Contours, 8 to 12 jtm, SPF2 Flowfield plus 2.5% Carbon, 
at t=5 sec 
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4.2.2 MNASA06/ASRM48-3 Inframetrics Thermal Image Camera Data 

The inframetrics 600 thermal image camera was used to obtain plume thermal images in 
the 8 to 12 /*m bandwidth for the ASRM48-3 test. These data were taken by Don Bryan (ED64) 
of NASA/MSFC and were digitized and provided to SECA, Inc. for this analysis (Ref. 4-6). 
The thermal image is the average of frames from 5.0 to 5.5 sec which provides a composite 
image at the 5.0 sec time frame. Plume temperature (°F) and radiance (watts/cm 2 /sr) contours 
from the Inframetrics 600 camera are shown in Figs. 4.12 and 4.13, respectively. The 
maximum measured radiance level of approximately 0.642 watts/cm 2 /sr shown in Fig. 4.13 is 
clearly at the plume centerline at the nozzle exit plane and the centerline radiance decreases 
approximately linearly for the first two nozzle exit radii downstream of the nozzle exit plane. 

The SIRRM-II predicted radiance contours for the ASRM48-3 test at t = 5 sec using a 
cycle 1.0 FDNS and cycle 2 SPF2 flowfields are shown in Figs. 4.14 and 4.15, respectively. 
Using the FDNS flowfield (Fig. 4.14), the predicted maximum plume centerline radiance of 
0.621 watts/cm 2 / sr occurs at approximately 5. 1 nozzle exit radii, while the predicted nozzle exit 
plane centerline radiance is 0.599 watts/cm 2 /sr which is 6.7% lower than the measurement 
maximum of 0.642 watts/cm 2 /sr at the same location. Using the SPF2 cycle 2 flowfield the 
maximum predicted radiance (Fig. 4.15) of 0.633 watts/cm 2 /sr occurs at the plume centerline 
at the nozzle exit plane and is 1.4% below the measured maximum value at the same location. 

A comparison of predicted and measured ASRM48-3 plume centerline radiance at t = 
5 sec is shown in Fig. 4.16. In this figure the SIRRM-II plume centerline radiance predictions 
using the FDNS, SPF2 cycle 2, and SPF2 cycle 1 flowfield methodologies are compared against 
the measured centerline radiance. The FDNSEL and SPF2 cycle 2 predictions are within 10% 
above the measurement from 0.5 to 7.5 nozzle exit radii downstream, and the SPF2 cycle 1 
prediction is within 13% below the measurement for the first four nozzle exit radii. All of the 
centerline radiance predictions shown in Fig. 4.16 fall within the generally accepted 
measurement accuracy of +20% for thermal imaging systems of this type. 
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Fig. 4.12 MNASA06/ASRM48-3 Inframetrics 600 Thermal Image Camera Temperature Contours (°F), 
8 to 12 /*m, Emissivity = 1.0, at t = 5.0 to 5.5 sec 
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MNASA06/ASRM48-3 SIRRM-II Predicted Radiance Contours, 8 to 12 /im, SPF2 Cycle 2 Flowfield, at t - 
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Fig. 4.16 MNASA06/ASRM48-3 Centerline Plume Radiance Comparison, 8 to 12 fim, at t = 5 sec 
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In the wavelength regions where MNASA thermal image data is available (2. 1 to 2.5 jtm, 
3.4 to 4.0 /xm, and 8 to 12 /*m), these data provide another valuable method of evaluating the 
existing solid motor plume radiation prediction capability. However, the bandwidth where solid 
particulate radiation dominates (1.0 to 2.0 jxm) has not been covered by the available thermal 
image measurements. Thermal image measurements in the 1.0 to 2.0 /*m bandwidth could 
provide a better understanding of the complete picture of solid motor radiation. 
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5.0 CONCLUSIONS AND RECOMMENDATIONS 


Based on the results of the research reported herein and on collateral work done by other 

investigators, the following conclusions and recommendations are offered. 

Conclusions: 

1. The Cycle 2 flowfield and particulate properties predictions when used with current 
optical property data for A1 2 0 3 provide acceptable estimates for SRM axisymmetric 
rocket plume base heating when used with a Monte Carlo radiation analysis. 

2. The generalized two-flux radiation model in the NOZZRAD code coupled with the view 
factor capability of the RAVFAC code provides a practical alternative to the Monte Carlo 
analysis and a convenient test vehicle for evaluating details of particulate/gas radiation 
analyses. 

3. The differences between the one-dimensional beam predictions made with the SIRRM and 
NOZZRAD codes requires further investigation. These differences are due to: (1) the 
method used for interpolation in A1 2 0 3 optical data tables, and (2) the method of coupling 
the particulate and gaseous radiation in specified wavelength regions. 

4. The method of spherical harmonics as developed herein appears to offer a practical 
alternative to the Monte Carlo analysis, but much more verification needs to be 
performed with this code before it is mature enough to be used for design purposes. 

5. The FDNS code with particle tracking has been shown to give comparable predictions 
to the two-phase RAMP code for SRM nozzle flows. Boundary conditions for the 
particle/gas mixture entering the nozzle are currently specified arbitrarily. A more 
rigorous analysis of the interior ballistics of the SRM grain should be developed and 
validated to improve the accuracy and utility of both of these flowfield codes. Additional 
analyses of SRM plumes need to be made with the FDNS code to validate this code to 
the same level as the SPF/2 code. Since the FDNS code conceptually treats 3- 
dimensional rocket plumes for which essentially no validation data exists, the 
FDNS/SPF/2 plume comparisons are currently the only method of validating the two- 
phase FDNS code as a plume code. 

6. The thermo-vision camera offers the potential for providing useful validation data for 
axisymmetric SRM plumes. However, the wavelength interval used for making the 
thermo-vision measurements should be carefully chosen so that the data can be accurately 
interpreted. For viewing the internal structure of the plumes, wavelengths which make 
the plume optically thin, which avoids regions of gas/particle interaction, and which 
avoids spectral intervals in which the particle optical properties undergo rapid changes 
should be chosen for the measuring system. Spectral intervals of 5-6 and of 7-8 /im 
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should be considered for internal viewing of plumes. To verify radiation analyses, 
spectral intervals of maximum heating should be observed (1-2 jtm). If the spectral 
intervals cannot be optimized because of instrumentation limitations, the thermo-vision 
camera data will be of limited value. 

7. The OD3P code provides mean particle size predictions that are consistent with the 
measurements taken by Sambamurthi during the MNASA and TEM test series. 

8. Preliminary exit plane radiation predictions for the MNASA motor indicate that the 
flowfield effects due to changes in combustion chamber geometry that occur during grain 
bumback can potentially explain the observed increase in measured radiation that occurs 
during a solid motor test firing. 

Recommendations: 

1. The A1 2 0 3 optical property data (N, and N 2 values) files for the Mie theory conversion 
to radiation properties (<r„ a b , and phase function parameters) should be optimized for 
radiation heating analysis. Gaseous/particle radiation interaction analysis should be 
further analyzed to include the recent improvements suggested by Reed (Ref. 5.1). 
These improvements should avoid the very narrow spectral interval and wide overall 
spectral region analyses required in the SIRRM code which were designed for plume 
signature analysis not base heating. 

2. The NOZZRAD code should be used as a preliminary design tool for radiation analysis. 
The REMCAR code for Monte Carlo predictions should be used where more detailed 
analyses are required. 

3. The method of spherical harmonics and extensions of the two-flux mode to make it a 
method of discrete coordinates should be further developed as alternative radiation 
analysis for future use. 

4. Additional two-phase FDNS code plume predictions should be compared to SPF/2 
predictions to obtain a validated CFD model for plume analyses. 

5. More interaction between plume analysts and thermo-vision instrumentation specialists 
should be made before new test programs are instituted. 
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